Questions about maniupating Sam format

[bbimber]

We are building an analysis pipeline for 454 sequence based genotyping. We compare each 454 read against a reference cDNA library and want to ask which sequences the read aligns against, and how many mismatches occur for each hit.

I have currently approached this by aligning the query reads against a ref library various aligners including bowtie and blat. The result is converted to SAM using perl, unless the aligner already outputs SAM.

From this SAM file, I am interested in extracting which reference sequences each read aligned against, and the SNPs associated with these alignments.

I have a few questions about manipulating SAM data:

1. From an alignment in SAM format, is there an established approach to annotate each alignment with the number of mismatches, their positions and/or identities (ie. G354C)? The normal CIGAR string does not distinguish match/mismatch although the doc says this is being considered.

2. Pileup does some of SNP detection discussed in point 1, but I did not see any information in the pileup report about which specifies query reads are associated with each SNP. If this information existed, it would be simple to use the result of pileup to annotate the SAM file. Can the read name(s) associated with each SNP be obtained from pileup or a similar tool?

3. This far, I have been aligning the experimental reads against a ref library of cDNA sequences. Tools like pileup will output reports in which each reference sequence is a 'chromosome' and list the reads that align to it. For my purposes, I am more interested in the reverse of this: which ref sequences align against each experimental read. Is there an established approach to transpose the ref/query sequences in a SAM file? It doesnt seem like it would be very hard in perl, but i thought i'd check to see if something already exists. Theoretically when I perform the alignment I could use my experimental data as the reference and align my cDNA library against it.

Thank you for any help.

[nilshomer]

Originally Posted by bbimber

We are building an analysis pipeline for 454 sequence based genotyping. We compare each 454 read against a reference cDNA library and want to ask which sequences the read aligns against, and how many mismatches occur for each hit.

I have currently approached this by aligning the query reads against a ref library various aligners including bowtie and blat. The result is converted to SAM using perl, unless the aligner already outputs SAM.

From this SAM file, I am interested in extracting which reference sequences each read aligned against, and the SNPs associated with these alignments.

I have a few questions about manipulating SAM data:

1. From an alignment in SAM format, is there an established approach to annotate each alignment with the number of mismatches, their positions and/or identities (ie. G354C)? The normal CIGAR string does not distinguish match/mismatch although the doc says this is being considered.

2. Pileup does some of SNP detection discussed in point 1, but I did not see any information in the pileup report about which specifies query reads are associated with each SNP. If this information existed, it would be simple to use the result of pileup to annotate the SAM file. Can the read name(s) associated with each SNP be obtained from pileup or a similar tool?

3. This far, I have been aligning the experimental reads against a ref library of cDNA sequences. Tools like pileup will output reports in which each reference sequence is a 'chromosome' and list the reads that align to it. For my purposes, I am more interested in the reverse of this: which ref sequences align against each experimental read. Is there an established approach to transpose the ref/query sequences in a SAM file? It doesnt seem like it would be very hard in perl, but i thought i'd check to see if something already exists. Theoretically when I perform the alignment I could use my experimental data as the reference and align my cDNA library against it.

Thank you for any help.

1.
The MD tag, if specified, can help you identify the mismatches without using the reference.

2. You can annotate SAM entries that span or include the SNP,.. The pileup does encode all the reads (see the last few columns), but they cannot be exactly identified by read name.

3. Depending on what tool you used, there will only be one mapping per read. Otherwise, others may be able to suggest tools.

[bbimber]

thanks for the reply. one more question:

an example pileup output, might look like:

seq1 272 T 24 ,.$.....,,.,.,...,,,.,..^+. <<<+;<<<<<<<<<<<=<;<;7<&

If I wanted to use this data to annotate a SAM file, I assume I simply find any lines where seq1 is the ref that span 272. What is the order of the characters ,.$.....,,.,.,...,,,.,..^+. relative to the lines of the SAM file? Will they match the order of whatever file was used to generate the pileup report (ie. the second character refers to the second line in the SAM where the ref sequence is seq2)? if true, it seems easy enough to find the associated read name and annotate them.

are there pre-existing tools or scripts that work to work backwards from a pileup output? is it common to work backwards from pileup output or does the pileup usually contain all the information people need?

[nilshomer]

Originally Posted by bbimber

thanks for the reply. one more question:

an example pileup output, might look like:

seq1 272 T 24 ,.$.....,,.,.,...,,,.,..^+. <<<+;<<<<<<<<<<<=<;<;7<&

If I wanted to use this data to annotate a SAM file, I assume I simply find any lines where seq1 is the ref that span 272. What is the order of the characters ,.$.....,,.,.,...,,,.,..^+. relative to the lines of the SAM file? Will they match the order of whatever file was used to generate the pileup report (ie. the second character refers to the second line in the SAM where the ref sequence is seq2)? if true, it seems easy enough to find the associated read name and annotate them.

are there pre-existing tools or scripts that work to work backwards from a pileup output? is it common to work backwards from pileup output or does the pileup usually contain all the information people need?

I don't know the order, you will have to consult the source code. Nonetheless, if two reads are aligned identically, you wont be able to identify them. It is not common to work backwards from the pileup alone, though going back into the SAM file is often necessary.