SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Tophat2: prepare unmapped.bam file for input into a tophat run on alternative genome danielsbrewer RNA Sequencing 10 05-20-2016 01:57 AM
Why reads in unmapped.bam still align to reference genome? SpreeFu Bioinformatics 7 09-28-2014 10:14 PM
How to extract unmapped reference sequence from BAM/SAM file.? Pinal Bioinformatics 2 03-06-2014 01:40 PM
Does accepted_hits.bam include unmapped.bam reads erichpeterson Bioinformatics 2 11-21-2012 08:08 AM
Obtaining reference identical reads from a BAM file Sakti Bioinformatics 2 05-17-2011 11:40 AM

Reply
 
Thread Tools
Old 09-16-2017, 12:43 PM   #1
Vca80553
Junior Member
 
Location: Stockholm

Join Date: Sep 2017
Posts: 5
Default Bam file with unmapped reads from another genome than reference

Hello everyone,

I started with bioinformatics 2 weeks ago, so maybe my question is a bit too easy for you, but I don't find any answer for it in other threads. I would appreciate a lot if you could help me out.

This is the case:

I mapped my paired end reads to my viral reference genome (Nextgenmap) and selected those that were mapped and paired (both pairs mapped only). Now I want to filter out, those reads that map also to human DNA (my contaminant). For that, I took the viral_mapped_paired_end_reads.bam, converted it to fastq (R1 and R2) and mapped it to hg19. In this bam file, I extracted the unmapped paired end. So, I assume that here I have the paired end reads that only mapped to virus before , but not to human now.
This bam file has the reads that I want, but unmapped to the human reference genome.

Now, how do I continue? I want to do coverage analysis for the viral genome for example. Can I use the unmapped bam file? or do I need to use the viral_mapped_paired_end_reads.bam and filter out the reads that mapped to human? If so, Is it done by extracting reads IDs? Or the IDs change depending on the reference genome?

Thanks a lot
Vca80553 is offline   Reply With Quote
Old 09-18-2017, 05:21 PM   #2
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,696
Default

I suggest you try BBSplit as in this thread. That will give you one fastq file for human reads, and one for viral reads. Then, map the viral output fastq file to the virus reference with a normal aligner (such as BBMap or Nextgenmap which you used previously). BBMap can output coverage data directly (well, BBSplit can too, actually...) using the covstats or basecov flags, if you want.

In answer to your last question, read IDs do not change with respect to reference genomes. However, unmapped bam files are not very useful and you can't use them for coverage analysis.
Brian Bushnell is offline   Reply With Quote
Old 09-19-2017, 05:27 AM   #3
Vca80553
Junior Member
 
Location: Stockholm

Join Date: Sep 2017
Posts: 5
Default

Dear Brian, Thanks a lot. I did as you suggested.
I also used BBMap (pileup.sh) for output coverage. I get the following results
Avg fold 53035,994
Length 7904
Ref_GC 0.0000
Covered % 100
Covered bases 7904
Plus reads 3190731
Minus reads 3190731
Read GC 0.395
Median_fold 22882
Std 24900.22

I assume I covered the whole genome. Do you happen to know whey the Ref_GC equals to 0.0000? I calculated GC content % of my reference genome and it is 36.51%. Thanks!
Vca80553 is offline   Reply With Quote
Old 09-19-2017, 01:16 PM   #4
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,696
Default

That's because you did not specify the reference file with "ref=" when you ran pileup.sh. It's not necessary; you only need it if you want the Ref_GC column to be correct.
Brian Bushnell is offline   Reply With Quote
Old 09-19-2017, 02:23 PM   #5
Vca80553
Junior Member
 
Location: Stockholm

Join Date: Sep 2017
Posts: 5
Default Ref_GC

Quote:
Originally Posted by Brian Bushnell View Post
That's because you did not specify the reference file with "ref=" when you ran pileup.sh. It's not necessary; you only need it if you want the Ref_GC column to be correct.
I wrote the following, but still didn't get it. Maybe something not right?

/home/sara/bbmap/pileup.sh in=1409_cat_sorted.bam "ref=/home/sara/HPV16Reference.fasta" out=1409.ref_stats1.txt hist=1409.ref1_histogram.txtll

Thanks
Vca80553 is offline   Reply With Quote
Old 09-19-2017, 03:48 PM   #6
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,696
Default

That's strange - I tested it and it works fine for me. No reference:

Code:
pileup.sh in=mapped.sam.gz stats=covstats.txt

cat covstats.txt

#ID     Avg_fold        Length  Ref_GC  Covered_percent Covered_bases   Plus_reads      Minus_reads     Read_GC Median_fold     Std_Dev
chr1    0.0908  249250621       0.0000  8.1241  20249466        76413   74517   0.4179  0       0.33
chr10   0.0960  135534747       0.0000  8.6332  11700973        43491   43228   0.4160  0       0.33
With reference:
Code:
pileup.sh in=mapped.sam.gz stats=covstats2.txt ref=hg19_main_mask_ribo_animal_allplant_allfungus.fa.gz

cat covstats2.txt

#ID     Avg_fold        Length  Ref_GC  Covered_percent Covered_bases   Plus_reads      Minus_reads     Read_GC Median_fold     Std_Dev
chr1    0.0908  249250621       0.4183  8.1241  20249466        76413   74517   0.4179  0       0.33
chr10   0.0960  135534747       0.4167  8.6332  11700973        43491   43228   0.4160  0       0.33
I tried various things including adding "hist=" and using "out=" instead of "stats=", and putting quotes around the reference flag, but was unable to replicate this. Are you sure you are looking at the correct output file, rather than an old one?
Brian Bushnell is offline   Reply With Quote
Reply

Tags
bam, other reference, unmapped reads

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 07:14 PM.


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