SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Truncated File: Trouble with SRA > Fastq > Sam DPCook Bioinformatics 0 12-10-2014 09:49 AM
smaller bam file given a bed file adrian Bioinformatics 2 09-17-2013 07:19 AM
Use SAM file to pull reads from FASTQ pbm13 Bioinformatics 5 06-29-2013 07:46 PM
Convert Sam file with incorrect fields to Fastq using HTSeq AnneBiton Bioinformatics 1 06-06-2013 07:03 AM
Split fastq into smaller files lorendarith Bioinformatics 10 12-13-2012 04:28 AM

Reply
 
Thread Tools
Old 09-30-2015, 08:59 AM   #1
scami
Member
 
Location: italy

Join Date: Sep 2010
Posts: 55
Default Sam file smaller than fastq

Hi Guys,

I hope it does not turn as a stupid question but this issue is driving me crazy. I have a reference genome A that contains chromosome I am not interested in. I used a python software to generate a new reference genome B with only some of the chromosomes of A. A and B have, as a way of example, Chr1 in common.
If I align all my reads on the complete reference genome A I get for Chr1 an average coverage of 10
If I align all my reads on chromosome B I get for Chr1 an average coverage of 2.
This does not make any sense! Moreover in the second case I obtain sam files that are smaller than the original reads fastq file, which I guess is also impossibile.
Any idea of such a strange result
I appreciate your help
scami is offline   Reply With Quote
Old 09-30-2015, 09:08 AM   #2
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,780
Default

How are you handling multimappers in both cases? Perhaps in case B reads are getting discarded if they multimap more than A. If all of your reads are not making it into the sam file (you are likely not including unmapped reads in your sam?) then it may explain the size result (though as rule of thumb don't worry about file sizes for NGS data as long as the reads are all accounted for).
GenoMax is online now   Reply With Quote
Old 09-30-2015, 09:14 AM   #3
scami
Member
 
Location: italy

Join Date: Sep 2010
Posts: 55
Default

Hi GenoMax

thanks for your reply. Parameters are exactly the same for the two alignments. I used bwa with the default setting except I used 0.05 as edit distance. I then used samtools to convert sam to bam, merge and sort bam file. In the past bwa produced a sam file that included all the reads of the fastq input file and I guess this has not changed. If so sam file, which contains both the sequence and the quality string, plus other fields should not be, as in my case, one third the length of the input fastq file.
Thanks.....
scami is offline   Reply With Quote
Old 09-30-2015, 09:20 AM   #4
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,780
Default

Have you checked your sam files using samtools idxstats?

See this recent thread as a reference for bam/sam file size issue (specially when the files are sorted).
GenoMax is online now   Reply With Quote
Old 09-30-2015, 10:42 PM   #5
scami
Member
 
Location: italy

Join Date: Sep 2010
Posts: 55
Default

Quote:
Originally Posted by GenoMax View Post
Have you checked your sam files using samtools idxstats?

See this recent thread as a reference for bam/sam file size issue (specially when the files are sorted).
Thanks for the advice. I ran the command you suggested and this is my result:
I divided the number of line of my fastq files in order to get the number of reads. I got a total of 187,595,025 reads. I calculated the sum of all mapped and unmapped reads in the bam file and obtained 65,536,000. I can not understand where all the other read went! Also in the idxstats output what is the difference between the ummapped reads for each chromosome and the unmapped reads at the end of the tab (the one with chromosome *):
This is my idxstats output:


chr1 23037639 2881935 72869
chr10 18140952 2435272 74191
chr11 19818926 2692338 74879
chr12 22702307 2797718 99784
chr13 24396255 2992318 96653
chr14 30274277 3751313 116119
chr15 20304914 2609991 112248
chr16 22053297 2630989 99628
chr17 17126926 2166178 70746
chr18 29360087 3792840 114974
chr19 24021853 3113695 120154
chr2 18779844 2507134 89904
chr3 19341862 2276796 82991
chr4 23867706 2918446 100922
chr5 25021643 3282028 116137
chr6 21508407 2523489 76380
chr7 21026613 2756721 79519
chr8 22385789 2801900 69991
chr9 23006712 3440921 151837
* 0 0 9344052

thanks again for your help!
scami is offline   Reply With Quote
Old 10-01-2015, 04:27 AM   #6
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,780
Default

Did you divide the number of lines by 4 to arrive at the ~187M number?

Here is a useful thread that explains the idxstats output: https://www.biostars.org/p/14569/
GenoMax is online now   Reply With Quote
Old 10-01-2015, 05:25 AM   #7
scami
Member
 
Location: italy

Join Date: Sep 2010
Posts: 55
Default

Quote:
Originally Posted by GenoMax View Post
Did you divide the number of lines by 4 to arrive at the ~187M number?

Here is a useful thread that explains the idxstats output: https://www.biostars.org/p/14569/
Yes I did. I also thought I messed up the reference genome file while removing the unwanted chromosomes with my python software. However afterwards I calculated the GC content of the chromosomes in the original and in the "stripped" reference file and they correspond exactly
scami 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 03:04 AM.


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