SEQanswers

Go Back   SEQanswers > Applications Forums > RNA Sequencing



Similar Threads
Thread Thread Starter Forum Replies Last Post
About Samtools view -L Brace Bioinformatics 2 03-29-2013 09:33 PM
Some confused from samtools mpileup. fabrice Bioinformatics 2 10-24-2012 06:02 PM
samtools view cmccabe Bioinformatics 0 07-21-2012 10:07 AM
samtools view seq_lover Bioinformatics 2 04-27-2012 12:22 PM
Samtools view michalkovac Bioinformatics 2 07-19-2011 06:25 AM

Reply
 
Thread Tools
Old 06-27-2014, 07:34 AM   #1
Rivalyn
Member
 
Location: Netherlands

Join Date: Apr 2014
Posts: 14
Default Confused about samtools view

This might be a silly question, but I'm confused about the output of samtools view. When I execute:

samtools view 1_BirA_mm10.bam | head

The output is this:

HWI-M01495:10:000000000-A5BGC:1:2107:12480:27069 337 1 30503670150M chrY 20418997 0 TCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTCCTTCTTCTTCTTCTTCTTCTTC HAGC/CHGFFHGHFHGCAGA0G/CG0FG/HAHHGCHFHGGHHFC?GFHGGFGEEAHHFCGGGCGCHFHHHEHGFFHFHFHHG1FHHGCHGFFFFHHHG3HHGAFHHEEGGCFADHGHG2HGAHFDGCABGCF4BA4FGFBFFFFCAAA@A AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:150YT:Z:UU NH:i:20 CC:Z:= CP:i:174010380 HI:i:3

As I understand it, this is not the Sam format. The Sam format should look like it does here:

http://samtools.github.io/hts-specs/SAMv1.pdf

Since you can use samtools view to convert from bam to sam, I'm trying to understand why the output is not in the sam format? Can someone explain the obvious thing I am apparently missing?
Rivalyn is offline   Reply With Quote
Old 06-27-2014, 07:49 AM   #2
WhatsOEver
Senior Member
 
Location: Germany

Join Date: Apr 2012
Posts: 215
Default

This is indeed a strange output. It looks like all fields are there but in the wrong order and with strange values. How did you produce this bam file / what mapper did you use? what version of samtools do you have?
WhatsOEver is offline   Reply With Quote
Old 06-27-2014, 07:52 AM   #3
Rivalyn
Member
 
Location: Netherlands

Join Date: Apr 2014
Posts: 14
Default

I used Tophat2, invoked like this:

tophat --num-threads 12 --library-type fr-unstranded -o ./tophat_trimmed_mm10/"${prefix%_*}" /home/kat/Data/Reference_sequences/Mus_musculus/UCSC/mm10/Sequence/Bowtie2Index/genome ./BIP_post_trim/Paired/Forward/"${b}" ./BIP_post_trim/Paired/Reverse/"${b2}"

Samtools version: 0.1.19-96b5f2294a
Rivalyn is offline   Reply With Quote
Old 06-27-2014, 08:10 AM   #4
WhatsOEver
Senior Member
 
Location: Germany

Join Date: Apr 2012
Posts: 215
Default

And there were no error printed at any point?
Have you build the bowtie2 index on your own (from the data structure, I would assume its an iGenomes download?!) - maybe it helps to re do that.
Do actually all entries look like that?
Does the sample dataset from tophat work for you?
WhatsOEver is offline   Reply With Quote
Old 06-27-2014, 08:21 AM   #5
Rivalyn
Member
 
Location: Netherlands

Join Date: Apr 2014
Posts: 14
Default

I can't see any errors in the terminal output. The terminal output is identical for all my biological replicates.

Yes, the reference file is from iGenomes. I didn't think it was necessary to build your own bowtie2 index if your reference is something as common as mouse.

Yes, all entries in the file look like that.

I ran the tophat test data quite some time ago. Let me see if I still have the output somewhere. Otherwise, I'll do it again.

4_Csde1_R1_paired.fq.gz and 4_Csde1_R2_paired.fq.gz match!

[2014-06-27 02:33:13] Beginning TopHat run (v2.0.9)
-----------------------------------------------
[2014-06-27 02:33:13] Checking for Bowtie
Bowtie version: 2.1.0.0
[2014-06-27 02:33:13] Checking for Samtools
Samtools version: 0.1.19.0
[2014-06-27 02:33:13] Checking for Bowtie index files (genome)..
[2014-06-27 02:33:13] Checking for reference FASTA file
[2014-06-27 02:33:13] Generating SAM header for /home/Data/Reference_sequences/Mus_musculus/UCSC/mm10/Sequence/Bowtie2Index/genome
format: fastq
quality scale: phred33 (default)
[2014-06-27 02:33:14] Preparing reads
left reads: min. length=50, max. length=150, 2067444 kept reads (43 discarded)
right reads: min. length=50, max. length=150, 2067265 kept reads (222 discarded)
[2014-06-27 02:34:48] Mapping left_kept_reads to genome genome with Bowtie2
[2014-06-27 02:40:24] Mapping left_kept_reads_seg1 to genome genome with Bowtie2 (1/6)
[2014-06-27 02:41:11] Mapping left_kept_reads_seg2 to genome genome with Bowtie2 (2/6)
[2014-06-27 02:42:03] Mapping left_kept_reads_seg3 to genome genome with Bowtie2 (3/6)
[2014-06-27 02:42:48] Mapping left_kept_reads_seg4 to genome genome with Bowtie2 (4/6)
[2014-06-27 02:43:29] Mapping left_kept_reads_seg5 to genome genome with Bowtie2 (5/6)
[2014-06-27 02:44:07] Mapping left_kept_reads_seg6 to genome genome with Bowtie2 (6/6)
[2014-06-27 02:44:33] Mapping right_kept_reads to genome genome with Bowtie2
[2014-06-27 02:50:00] Mapping right_kept_reads_seg1 to genome genome with Bowtie2 (1/6)
[2014-06-27 02:50:45] Mapping right_kept_reads_seg2 to genome genome with Bowtie2 (2/6)
[2014-06-27 02:51:42] Mapping right_kept_reads_seg3 to genome genome with Bowtie2 (3/6)
[2014-06-27 02:52:31] Mapping right_kept_reads_seg4 to genome genome with Bowtie2 (4/6)
[2014-06-27 02:53:13] Mapping right_kept_reads_seg5 to genome genome with Bowtie2 (5/6)
[2014-06-27 02:53:49] Mapping right_kept_reads_seg6 to genome genome with Bowtie2 (6/6)
[2014-06-27 02:54:12] Searching for junctions via segment mapping
[2014-06-27 02:57:15] Retrieving sequences for splices
[2014-06-27 02:58:33] Indexing splices
[2014-06-27 02:58:48] Mapping left_kept_reads_seg1 to genome segment_juncs with Bowtie2 (1/6)
[2014-06-27 02:59:03] Mapping left_kept_reads_seg2 to genome segment_juncs with Bowtie2 (2/6)
[2014-06-27 02:59:20] Mapping left_kept_reads_seg3 to genome segment_juncs with Bowtie2 (3/6)
[2014-06-27 02:59:36] Mapping left_kept_reads_seg4 to genome segment_juncs with Bowtie2 (4/6)
[2014-06-27 02:59:52] Mapping left_kept_reads_seg5 to genome segment_juncs with Bowtie2 (5/6)
[2014-06-27 03:00:08] Mapping left_kept_reads_seg6 to genome segment_juncs with Bowtie2 (6/6)
[2014-06-27 03:00:22] Joining segment hits
[2014-06-27 03:02:20] Mapping right_kept_reads_seg1 to genome segment_juncs with Bowtie2 (1/6)
[2014-06-27 03:02:36] Mapping right_kept_reads_seg2 to genome segment_juncs with Bowtie2 (2/6)
[2014-06-27 03:02:52] Mapping right_kept_reads_seg3 to genome segment_juncs with Bowtie2 (3/6)
[2014-06-27 03:03:08] Mapping right_kept_reads_seg4 to genome segment_juncs with Bowtie2 (4/6)
[2014-06-27 03:03:23] Mapping right_kept_reads_seg5 to genome segment_juncs with Bowtie2 (5/6)
[2014-06-27 03:03:38] Mapping right_kept_reads_seg6 to genome segment_juncs with Bowtie2 (6/6)
[2014-06-27 03:03:52] Joining segment hits
[2014-06-27 03:05:52] Reporting output tracks
-----------------------------------------------
[2014-06-27 02:33:13] A summary of the alignment counts can be found in ./tophat_trimmed_mm10/4_Csde1/align_summary.txt
[2014-06-27 02:33:13] Run complete: 00:49:30 elapsed

Last edited by Rivalyn; 06-27-2014 at 09:08 AM.
Rivalyn is offline   Reply With Quote
Old 06-27-2014, 08:27 AM   #6
Rivalyn
Member
 
Location: Netherlands

Join Date: Apr 2014
Posts: 14
Default

I found the output from the tophat test data. All of the reads look like this:

test_mRNA_3_187_51 99 test_chromosome 53 50 75M = 163 185 TACTATTTGACTAGACTGGAGGCGCTTGCGACTGAGCTAGGACGTGCCACTACGGGGATGACGACTCGGACTACG IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII AS:i:-12 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:6C59A8 YT:Z:UU NH:i:1
Rivalyn is offline   Reply With Quote
Old 07-01-2014, 04:49 AM   #7
WhatsOEver
Senior Member
 
Location: Germany

Join Date: Apr 2012
Posts: 215
Default

Well, the test data and your run-log look perfectly fine...

However, you're versions of tophat and bowtie are quite outdated - the current releases are 2.0.12 and 2.2.3. It doesn't make much sense to me how this could have resulted in an output like you showed, but it's worth trying to update them...

I always build the indices myself as only then you can be sure that no strange version issues occur.

Btw: Have you tried a different mapper?
WhatsOEver is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 10:25 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO