Note
Make sure you’re running in screen!
Start with the QC’ed files from 1. Quality Trimming and Filtering Your Sequences or copy them into a working directory;
Normalize everything to a coverage of 20, starting with the (more valuable) PE reads; keep pairs using ‘-p’
cd /mnt/work
normalize-by-median.py -p -k 20 -C 20 -N 4 -x 1e9 --savetable normC20k20.kh *.pe.qc.fq.gz
and continuing into the (less valuable but maybe still useful) SE reads
normalize-by-median.py -C 20 --loadtable normC20k20.kh --savetable normC20k20.kh *.se.qc.fq.gz
This produces a set of ‘.keep’ files, as well as a normC20k20.kh database file.
Use ‘filter-abund’ to trim off any k-mers that are abundance-1 in high-coverage reads. The -V option is used to make this work better for variable coverage data sets:
filter-abund.py -V normC20k20.kh *.keep
This produces .abundfilt files containing the trimmed sequences.
The process of error trimming could have orphaned reads, so split the PE file into still-interleaved and non-interleaved reads
for i in *.pe.qc.fq.gz.keep.abundfilt
do
extract-paired-reads.py $i
done
This leaves you with PE files (.pe.qc.fq.gz.keep.abundfilt.pe) and two sets of SE files (.se.qc.fq.gz.keep.abundfilt and .pe.qc.fq.gz.keep.abundfilt.se). (Yes, the naming scheme does make sense. Trust me.)
You can now remove the normC20k20.kh file, too
rm normC20k20.kh
Now that we’ve eliminated many more erroneous k-mers, let’s ditch some more high-coverage data. First, normalize the paired-end reads
normalize-by-median.py -C 5 -k 20 -N 4 -x 1e9 --savetable normC5k20.kh -p *.pe.qc.fq.gz.keep.abundfilt.pe
and then do the remaining single-ended reads
normalize-by-median.py -C 5 --savetable normC5k20.kh --loadtable normC5k20.kh *.pe.qc.fq.gz.keep.abundfilt.se *.se.qc.fq.gz.keep.abundfilt
Now let’s tidy things up. Here are the paired files (kak = keep/abundfilt/keep)
for i in *.abundfilt.pe
do
newfile=$(basename $i .pe.qc.fq.gz.keep.abundfilt.pe)
echo newfile is $newfile
gzip -c $i > $newfile.pe.kak.qc.fq.gz
done
and for the single-ended files
for i in *.se.qc.fq.gz.keep.abundfilt
do
pe_orphans=$(basename $i .se.qc.fq.gz.keep.abundfilt).pe.qc.fq.gz.keep.abundfilt.se
newfile=$(basename $i .se.qc.fq.gz.keep.abundfilt).se.kak.qc.fq.gz
cat $i $pe_orphans | gzip -c > $newfile
done
You can now remove all of these various files:
*.pe.qc.fq.gz.keep
*.pe.qc.fq.gz.keep.abundfilt
*.pe.qc.fq.gz.keep.abundfilt.pe
*.pe.qc.fq.gz.keep.abundfilt.pe.keep
*.pe.qc.fq.gz.keep.abundfilt.se
*.pe.qc.fq.gz.keep.abundfilt.se.keep
by typing
rm *.keep *.abundfilt *.pe *.se
If you are not doing partitioning (see 3. Partitioning), you may also want to remove the k-mer hash tables:
rm *.kh
If you are running partitioning, you can remove the normC20k20.kh file:
rm normC20k20.kh
but you will need the normC5k20.kh file.
Try running
/usr/local/share/khmer/sandbox/readstats.py *.kak.qc.fq.gz *.?e.qc.fq.gz
after a long wait, you’ll see:
---------------
861769600 bp / 8617696 seqs; 100.0 average length -- SRR606249.pe.qc.fq.gz
79586148 bp / 802158 seqs; 99.2 average length -- SRR606249.se.qc.fq.gz
531691400 bp / 5316914 seqs; 100.0 average length -- SRR606249.pe.qc.fq.gz
89903689 bp / 904157 seqs; 99.4 average length -- SRR606249.se.qc.fq.gz
173748898 bp / 1830478 seqs; 94.9 average length -- SRR606249.pe.kak.qc.fq.gz
8825611 bp / 92997 seqs; 94.9 average length -- SRR606249.se.kak.qc.fq.gz
52345833 bp / 550900 seqs; 95.0 average length -- SRR606249.pe.kak.qc.fq.gz
10280721 bp / 105478 seqs; 97.5 average length -- SRR606249.se.kak.qc.fq.gz
---------------
This shows you how many sequences were in the original QC files, and how many are left in the ‘kak’ files. Not bad – considerably more than 80% of the reads were eliminated in the kak!
Next: 3. Partitioning