[nilshomer]

Hello,

I was wondering if there was a score associated with each read in the SAM file, that would give an indication on the strength of the match between the read and the subject sequence. The CIGAR string helps to some extent, but since "M" denotes match or mismatch, I was wondering if there was a way to differentiate between the two. Sort of like an E-value or a blast score.

Abhijit

Probably a good idea to create a new thread (this one is getting long!).
See the mapping quality field.

Hi guys,

Had a quick question regarding the SAM- CIGAR column. I understand the M attribute designates both matches and mismatches. Is there a way to get at the literal number of mismatches without resorting to comparing the tags or sequence to the reference using the SAM format? Sorry if this question is a re-post. I tried searching, but couldn't find anything.

Probably a good idea to create a new thread (this one is getting long!).
Try the NM optional tag if it is available (aligner specific).

[JohnK]

Hi guys,

Had a quick question regarding the SAM- CIGAR column. I understand the M attribute designates both matches and mismatches. Is there a way to get at the literal number of mismatches without resorting to comparing the tags or sequence to the reference using the SAM format? Sorry if this question is a re-post. I tried searching, but couldn't find anything.

[gen2prot]

Hello,

I was wondering if there was a score associated with each read in the SAM file, that would give an indication on the strength of the match between the read and the subject sequence. The CIGAR string helps to some extent, but since "M" denotes match or mismatch, I was wondering if there was a way to differentiate between the two. Sort of like an E-value or a blast score.

Abhijit

[gen2prot]

Hello nilshomer,

I specified the SQ field but I still get an error, probably because the reference sequences are drosophila gene sequences. Therefore the names of the reference sequences are different unless 2 or more reads match to the same gene. I cannot convert the names to a single reference sequence name, since I will loose information. I am stuck. Maybe I need to do the traditional perl sort (very time consuming for a 6GB file). Any better way of doing this?

Thanks
Abhijit

[nilshomer]

Hello,

I gave the following command using the Picard tool Sortsam

java -jar SortSam.jar I=../test/testsam.sam O=../test/sortedtest.sam SO=queryname

My Input file looks like this:

@HD VN:1.0 SO:sorted
@PG ID:TopHat VN:1.0.13 CL:/share/apps/bin/tophat -o ./s1 --solexa1.3-quals -p 2 GeneIndex /home/asanyal/data/Flydata/Exp_100423/100423_HWI-EAS313_0001_61G2CAAXX.birchlerj/s_1_sequence.txt
HWI-EAS313_0001:1:80:8942:6680#0 0 FBgn0000003 1 3 42M * 0 0 CGGACTGGAAGGTTGGCAGCTTCTGTAATCACGCTTCTGTGA GGGFGGGGGGGFFEGGGGGGGFGGGGGGDDGFGGGGGGGGFE NM:i:2
HWI-EAS313_0001:1:108:8254:11808#0 0 FBgn0000003 9 3 42M * 0 0 AAGGTTGGCAGCTTCTGTAATCACGCTTCTGTGAGGTCTGAT C::?>ACCCCD?EDEB=EEEEEECEE?:E??@C@CEBED=4? NM:i:0

However, I get the following error message.

Exception in thread "main" net.sf.samtools.SAMFormatException: Error parsing text SAM file. Empty sequence dictionary.; Line 3
Line: HWI-EAS313_0001:1:80:8942:6680#0 0 FBgn0000003 1 3 42M * 0 0 CGGACTGGAAGGTTGGCAGCTTCTGTAATCACGCTTCTGTGA GGGFGGGGGGGFFEGGGGGGGFGGGGGGDDGFGGGGGGGGFE NM:i:2

Do I have to give the program the reference sequences? Or do I need to create a sequence dictionary using CreateSequenceDictionary

Thanks
Abhijit

There are no "SQ" fields in your SAM file. You could try giving it the reference sequence if they are not present.

[gen2prot]

Hello,

I gave the following command using the Picard tool Sortsam

java -jar SortSam.jar I=../test/testsam.sam O=../test/sortedtest.sam SO=queryname

My Input file looks like this:

@HD VN:1.0 SO:sorted
@PG ID:TopHat VN:1.0.13 CL:/share/apps/bin/tophat -o ./s1 --solexa1.3-quals -p 2 GeneIndex /home/asanyal/data/Flydata/Exp_100423/100423_HWI-EAS313_0001_61G2CAAXX.birchlerj/s_1_sequence.txt
HWI-EAS313_0001:1:80:8942:6680#0 0 FBgn0000003 1 3 42M * 0 0 CGGACTGGAAGGTTGGCAGCTTCTGTAATCACGCTTCTGTGA GGGFGGGGGGGFFEGGGGGGGFGGGGGGDDGFGGGGGGGGFE NM:i:2
HWI-EAS313_0001:1:108:8254:11808#0 0 FBgn0000003 9 3 42M * 0 0 AAGGTTGGCAGCTTCTGTAATCACGCTTCTGTGAGGTCTGAT C::?>ACCCCD?EDEB=EEEEEECEE?:E??@C@CEBED=4? NM:i:0

However, I get the following error message.

Exception in thread "main" net.sf.samtools.SAMFormatException: Error parsing text SAM file. Empty sequence dictionary.; Line 3
Line: HWI-EAS313_0001:1:80:8942:6680#0 0 FBgn0000003 1 3 42M * 0 0 CGGACTGGAAGGTTGGCAGCTTCTGTAATCACGCTTCTGTGA GGGFGGGGGGGFFEGGGGGGGFGGGGGGDDGFGGGGGGGGFE NM:i:2

Do I have to give the program the reference sequences? Or do I need to create a sequence dictionary using CreateSequenceDictionary

Thanks
Abhijit

[gen2prot]

Thank you.

[nilshomer]

Hello nilshomer,

I downloaded picard. I have the .jar files on MAC osx 10.6. Yet these jar files won't open. I have them saved on the Desktop. How do I run it?

Thanks
Abhijit

I am assuming you have familiarity with the Terminal and a Unix-based environment. If this is wrong, you need to become familiar with these environments (search this site for recommended books and tutorials). I cannot teach you how to use the Terminal and such basic questions.

Use the command for the respective jar:

java -jar SortSam.jar

[gen2prot]

Hello nilshomer,

I downloaded picard. I have the .jar files on MAC osx 10.6. Yet these jar files won't open. I have them saved on the Desktop. How do I run it?

Thanks
Abhijit

[nilshomer]

Hello All,

Does anybody know how can I sort the .sam file on the basis of the first column? That is the column containing the unique read identifiers? Right now its sorted on the 3rd.

Thanks
Abhijit

SAMtools and Picard will both sort by read name. See their documentation.

[gen2prot]

Hello All,

Does anybody know how can I sort the .sam file on the basis of the first column? That is the column containing the unique read identifiers? Right now its sorted on the 3rd.

Thanks
Abhijit

[bair]

denovo on sam format

Dear all,

I have alignment results in bam file which includes pair-end, mate-pair reads in different length (101 and 35, 36bp). Does anybody know that Soap or other denovo program can handle with bam format directly or I have to use the raw reads files?

Many thanks!

[menenuh]

non-unique reads

Hello,
In my sam file I have both unique and non-unique reads. What happens to non-unique reads when I call SNPs from the sam file? Are they included in the SNP calling process?

thanks

[RockChalkJayhawk]

Perfect!

Thanks Nils! Youre the best!

[nilshomer]

So in that case, my MRNM does not equal "=" OR MRNM equals "=" and the difference between POS and MPOS > 1 million.

Is this correct?

Perfect!