Getting reads from lastal results?

[kga1978]

Hi All,

I have a bunch of paired-end reads and I aligned them against a custom database (viruses) using LASTAL:

lastal -Q1 /idi/sabetilab/kandersen/references/blast/arena G771_2.reads1.sub.fastq | maf-sort.sh -n2 > G771_2.reads1.sub.maf
lastal -Q1 /idi/sabetilab/kandersen/references/blast/arena G771_2.reads2.sub.fastq | maf-sort.sh -n2 > G771_2.reads2.sub.maf

I then used the "last-pair-probs" script to create a single maf file containing the entries:

last-pair-probs.py G771_2.reads1.sub.maf G771_2.reads2.sub.maf > G771_2.sub.maf

Followed by conversion into SAM format (to get the reads):

maf-convert.py sam G771_2.sub.maf > G771_2.sub.sam

So far, so good, except I think the SAM file doesn't correctly assign the two mates. Here's what the resultant SAM file looks like:

@HD	VN:1.3	SO:unknown
D0N2CACXX120229:2:1101:21311:18913/1	0	LASV-G771-S-Sierra_Leone-2010H	14	24	101M	*	0	0	CCTAGGCATTTTTGGTTGCGCAATTCAAGTGT
CCTATTTAAAATGGGACAGATAGTGACATTCTTCCAGGAAGTGCCTCATGTAATAGAAGAGGTGATGAA	@@@FFFDEHHHHFEGIJGBFHIIGG>DD>GIGGIDDFEEF@>DG>DC3BF;=<@FC@F@FGGC;D@=>AEA;?)7?.;@3>CE;>>;(
353559@A#####	NM:i:0	AS:i:584
D0N2CACXX120229:2:1101:21311:18913/2	16	LASV-G771-S-Sierra_Leone-2010H	646	24	101M	*	0	0	TGGTATTTACATTGCTCTTGACTCAGGCCGTG
ACCGGTGGGACTGTATTATGACTAGTTATCAATATCTGATAATCCAAAATACGACCTGGGAAGATCACT	################@;@>5==A?DFEECA/@E@BF@=B0?4B4DFDB<??>HCFB<<D<>G?9D?ED?GFEE@GEFJIEEGGBAHF

As you can see, the two mates are in there, but when I run picard to convert to FASTQ format (required for downstream analyses), then I only get the first mate ("G771_2.sub.sam.reads1.fastq " containing all /1 and /2 reads), but not the second ("G771_2.sub.sam.reads.reads2.fastq" - this file is empty):

SamToFastq.jar INPUT=G771_2.sub.sam FASTQ=G771_2.sub.sam.reads1.fastq SECOND_END_FASTQ=G771_2.sub.sam.reads.reads2.fastq VALIDATION_STRINGENCY=SILENT

@D0N2CACXX120229:2:1101:21311:18913/1
CCTAGGCATTTTTGGTTGCGCAATTCAAGTGTCCTATTTAAAATGGGACAGATAGTGACATTCTTCCAGGAAGTGCCTCATGTAATAGAAGAGGTGATGAA
+
@@@FFFDEHHHHFEGIJGBFHIIGG>DD>GIGGIDDFEEF@>DG>DC3BF;=<@FC@F@FGGC;D@=>AEA;?)7?.;@3>CE;>>;(353559@A#####
@D0N2CACXX120229:2:1101:21311:18913/2
AGTGATCTTCCCAGGTCGTATTTTGGATTATCAGATATTGATAACTAGTCATAATACAGTCCCACCGGTCACGGCCTGAGTCAAGAGCAATGTAAATACCA
+
@?@DDDDDHDCDFFHABGGEEIJFEG@EEFG?DE?D9?G><D<<BFCH>??<BDFD4B4?0B=@FB@E@/ACEEFD?A==5>@;@################

Any thoughts on how to resolve this? In other word, how to take a bunch of paired-end reads, align them against a database using LASTAL followed by extraction of the aligned (or for other purposes - unaligned) reads in the correct mate pairs?

Thanks very much in advance.

[maubp]

...
So far, so good, except I think the SAM file doesn't correctly assign the two mates. Here's what the resultant SAM file looks like:

@HD	VN:1.3	SO:unknown
D0N2CACXX120229:2:1101:21311:18913/1	0	LASV-G771-S-Sierra_Leone-2010H	14	24	101M	*	0	0	CCTAGGCATTTTTGGTTGCGCAATTCAAGTGT
CCTATTTAAAATGGGACAGATAGTGACATTCTTCCAGGAAGTGCCTCATGTAATAGAAGAGGTGATGAA	@@@FFFDEHHHHFEGIJGBFHIIGG>DD>GIGGIDDFEEF@>DG>DC3BF;=<@FC@F@FGGC;D@=>AEA;?)7?.;@3>CE;>>;(
353559@A#####	NM:i:0	AS:i:584
D0N2CACXX120229:2:1101:21311:18913/2	16	LASV-G771-S-Sierra_Leone-2010H	646	24	101M	*	0	0	TGGTATTTACATTGCTCTTGACTCAGGCCGTG
ACCGGTGGGACTGTATTATGACTAGTTATCAATATCTGATAATCCAAAATACGACCTGGGAAGATCACT	################@;@>5==A?DFEECA/@E@BF@=B0?4B4DFDB<??>HCFB<<D<>G?9D?ED?GFEE@GEFJIEEGGBAHF

Yes, that is not a proper SAM file if they are meant to be paired reads. The QNAME (column 1) should be just D0N2CACXX120229:2:1101:21311:18913 for these entries, and the first/second part encoded in the FLAG (column 2).

Something like this:

@HD	VN:1.3	SO:unknown
D0N2CACXX120229:2:1101:21311:18913	99	LASV-G771-S-Sierra_Leone-2010H	14	24	101M	*	0	0	CCTAGGCATTTTTGGTTGCGCAATTCAAGTGT
CCTATTTAAAATGGGACAGATAGTGACATTCTTCCAGGAAGTGCCTCATGTAATAGAAGAGGTGATGAA	@@@FFFDEHHHHFEGIJGBFHIIGG>DD>GIGGIDDFEEF@>DG>DC3BF;=<@FC@F@FGGC;D@=>AEA;?)7?.;@3>CE;>>;(
353559@A#####	NM:i:0	AS:i:584
D0N2CACXX120229:2:1101:21311:18913	147	LASV-G771-S-Sierra_Leone-2010H	646	24	101M	*	0	0	TGGTATTTACATTGCTCTTGACTCAGGCCGTG
ACCGGTGGGACTGTATTATGACTAGTTATCAATATCTGATAATCCAAAATACGACCTGGGAAGATCACT	################@;@>5==A?DFEECA/@E@BF@=B0?4B4DFDB<??>HCFB<<D<>G?9D?ED?GFEE@GEFJIEEGGBAHF

This is how I calculated the FLAG values (in Python), assuming they are properly paired and only the second is mapped on the reverse strand:

>>> 0x1 + 0x2 + 0x40 + 0x20
99
>>> 0x1 + 0x2 + 0x80 + 0x10
147

[kga1978]

I suspected it would be something like this - crumbs. I'll see if I can get around the problem some other way. Thanks maubp