Filter out low coverage contigs
0
0
Entering edit mode
2.0 years ago

I ran the Spades software (https://github.com/ablab/spades), which outputs contigs.fasta

for f in `ls -1 *_1.fq | sed 's/_1.fq//’`;
    do spades.py -o . --mismatch-correction -1 $f\_1.fq -2 $f\_2.fq -t 20;
done

Within the same directory, I created and executed a bash script filter_contigs.sh:

for file in contigs.fasta; do
        grep -F ">" $file | sed -e 's/_/ /g' |sort -nrk 6 |awk '$6>=2.0 && $4>=500 {print $0}'|
        sed -e 's/ /_/g'|sed -e 's/>//g'>$file.txt
        echo sequences to keep
        wc -l $file.txt
        echo running fastagrep.pl
        ../fastagrep.pl -f $file.txt $file > HCov.$file
        echo sequences kept
        grep -c ">" HCov.$file
done

When I executed filter_contigs.sh, it fails to capture/keep any contigs. Why?

sequences to keep
0 contigs.fasta.txt
running fastagrep.pl
sequences kept
0
bash contigs • 1.4k views
ADD COMMENT
0
Entering edit mode

what is the output of

grep -F ">" contigs.fasta

?

ADD REPLY
0
Entering edit mode

hmmm I didn't realize contigs.fasta is empty. Perhaps I will rerun spades. I think --mismatch-correction option is not available for my Spades' version.

ADD REPLY
0
Entering edit mode

Maybe --only-error-correction ?

ADD REPLY
0
Entering edit mode

I want to both correct and assemble the reads (I removed the previous error correction steps).

ADD REPLY
0
Entering edit mode

There is no switch needed for both error-correction and assembly, as that is a default behavior.

ADD REPLY
0
Entering edit mode

But if I use --only-error-correction, doesn't it mean that it won't perform the assembly? If I run spades.py -o . -1 read_1 -2 read_2 -t 20, it would perform both error correction and assembly, right?

ADD REPLY
0
Entering edit mode

A word of caution; it's been a while since I've used SPAdes so maybe this has changed but I'm pretty sure the "coverage" value in the header of SPAdes contigs is not sequence coverage, it is kmer coverage.

ADD REPLY

Login before adding your answer.

Traffic: 2018 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6