SEQanswers

Go Back   SEQanswers > General



Similar Threads
Thread Thread Starter Forum Replies Last Post
the low percentage of mapped reads wangli RNA Sequencing 2 04-04-2012 05:56 AM
Calculating Percentage of n's for Each Transcript in FASTA File KDS Bioinformatics 7 06-07-2011 11:06 PM
a question for the percentage of aligned reads chenwb SOLiD 1 03-18-2011 07:58 AM
Percentage of mapped reads ? zack80.liu Bioinformatics 6 03-01-2011 09:08 AM
low percentage of reads mapped rahilsethi SOLiD 3 09-13-2010 06:01 AM

Reply
 
Thread Tools
Old 04-08-2015, 01:01 AM   #1
kaps
Member
 
Location: Uganda

Join Date: Jan 2015
Posts: 71
Default Calculating percentage of reads aligning to a subject

I have done a blast on virus contigs and obtained hits matching to virus (strains/isolates) in the database. I would like to calculate the percentage of reads that are aligning to the viruses in the database. Can someone guide me on how to do this?
kaps is offline   Reply With Quote
Old 04-08-2015, 02:37 AM   #2
Michael.Ante
Senior Member
 
Location: Vienna

Join Date: Oct 2011
Posts: 121
Default

You can download the virus sequences in fasta format and use e.g. Bowtie2 to align the reads locally. Therefore, you need to build an index first. The log output of Bowtie2 tells you haw many reads mapped.
After aligning the reads, you can use samtools to get some statistics (e.g. samtools idxstats).
Michael.Ante is offline   Reply With Quote
Old 04-11-2015, 02:38 AM   #3
kaps
Member
 
Location: Uganda

Join Date: Jan 2015
Posts: 71
Default

Thanks Michael,

I have not used Bowtie/ samtools before. How do I start off?
kaps is offline   Reply With Quote
Old 04-11-2015, 03:39 AM   #4
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,975
Default

1. What format are your blast results in (html, xml, text)? You may be able to parse that result file if all you want to know is how many sequences hit a "virus".

2. If you did the blast locally do you have a sequence file with all "virus" sequences available? You will be able to use that file as an input for bowtie2 and follow the path @Michael.Ante suggested.

3. Are you comfortable using command line (e.g. linux) applications?
GenoMax is offline   Reply With Quote
Old 04-11-2015, 04:57 AM   #5
kaps
Member
 
Location: Uganda

Join Date: Jan 2015
Posts: 71
Default

1. blast results are in txt format
2. yes. sequence file is available ( database file?)
3. am fairly comfortable with command line
4. I would prefer the command prompt option as opposed to logging on the cluster (my internet is erratic)

Last edited by kaps; 04-11-2015 at 05:05 AM.
kaps is offline   Reply With Quote
Old 04-14-2015, 12:23 AM   #6
Michael.Ante
Senior Member
 
Location: Vienna

Join Date: Oct 2011
Posts: 121
Default

Quote:
Originally Posted by kaps View Post
Thanks Michael,

I have not used Bowtie/ samtools before. How do I start off?
Hi Kaps,

you should have a look at the Bowtie2 homepage. There, it is explained in detail how the programs work. At the end of the manual is a "Lambda phage example", which has quite an overlap to your problem. It also has a SAMtools downstream section...

Cheers,
Michael
Michael.Ante is offline   Reply With Quote
Old 04-19-2015, 11:57 PM   #7
kaps
Member
 
Location: Uganda

Join Date: Jan 2015
Posts: 71
Default

Hi, Michael

Thanks for pointing this to me!
kaps is offline   Reply With Quote
Old 04-22-2015, 01:54 AM   #8
kaps
Member
 
Location: Uganda

Join Date: Jan 2015
Posts: 71
Default

Quote:
Originally Posted by Michael.Ante View Post
You can download the virus sequences in fasta format and use e.g. Bowtie2 to align the reads locally. Therefore, you need to build an index first. The log output of Bowtie2 tells you haw many reads mapped.
After aligning the reads, you can use samtools to get some statistics (e.g. samtools idxstats).
Hello Michael, when I try samtools idxstats,
I am getting a comment as below;

samtools idxstats lib4seq.sorted.bam
[bam_idxstats] fail to load the index.

what could be the problem?
kaps is offline   Reply With Quote
Old 04-22-2015, 02:42 AM   #9
Michael.Ante
Senior Member
 
Location: Vienna

Join Date: Oct 2011
Posts: 121
Default

Have you created the index with
Code:
samtools index lib4seq.sorted.bam
?
If yes how does your bam-dile header looks like?
Code:
 samtools view -H lib4seq.sorted.bam
Michael.Ante is offline   Reply With Quote
Old 04-22-2015, 06:05 AM   #10
kaps
Member
 
Location: Uganda

Join Date: Jan 2015
Posts: 71
Default

I had not created the index,
I can now see the statistics in the index file!

Thanks
kaps is offline   Reply With Quote
Old 04-23-2015, 05:43 AM   #11
kaps
Member
 
Location: Uganda

Join Date: Jan 2015
Posts: 71
Default

Quote:
Originally Posted by Michael.Ante View Post
You can download the virus sequences in fasta format and use e.g. Bowtie2 to align the reads locally. Therefore, you need to build an index first. The log output of Bowtie2 tells you haw many reads mapped.
After aligning the reads, you can use samtools to get some statistics (e.g. samtools idxstats).
In a case where my index file has several sequences for different strains/isolates of the same virus which may be treated as duplicates, how do I restrict bowtie2 to do the alignment once?
kaps is offline   Reply With Quote
Old 05-11-2015, 11:04 PM   #12
kaps
Member
 
Location: Uganda

Join Date: Jan 2015
Posts: 71
Default

Quote:
Originally Posted by Michael.Ante View Post
You can download the virus sequences in fasta format and use e.g. Bowtie2 to align the reads locally. Therefore, you need to build an index first. The log output of Bowtie2 tells you haw many reads mapped.
After aligning the reads, you can use samtools to get some statistics (e.g. samtools idxstats).
Hello,

After getting the samtools idxstats (on number of mapped vs unmapped reads), is it possible to extract/select reads that mapped from the raw read files/query? how is it done?
kaps is offline   Reply With Quote
Old 05-12-2015, 03:20 AM   #13
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,975
Default

If you had used the "--un-conc and --al-conc" options (http://bowtie-bio.sourceforge.net/bo...output-options) the unmapped reads could have been written to separate files when you did the alignment.


1. You could repeat bowtie2 alignment with above parameters added to your original list (easier) OR
2. Identify read ID's of sequences that mapped and use a tool like seqtk to extract the mapped reads (e.g. seqtk subseq in.fq name.lst > out.fq)


Use @Michael.Ante's easy suggestion below

Last edited by GenoMax; 05-12-2015 at 04:34 AM.
GenoMax is offline   Reply With Quote
Old 05-12-2015, 04:24 AM   #14
Michael.Ante
Senior Member
 
Location: Vienna

Join Date: Oct 2011
Posts: 121
Default

You can use samtools view to extract the mapped/unmapped reads by filtering the 'unmapped' flag:

Code:
samtools view -F 4 -bh lib4seq.sorted.bam > lib4seq.sorted.mapped.bam
samtools view -f 4 -bh lib4seq.sorted.bam > lib4seq.sorted.unmapped.bam
Samtools view will help a lot; just have a look at some tutorials, for instance Dave's wiki
Michael.Ante is offline   Reply With Quote
Old 05-19-2015, 02:38 AM   #15
kaps
Member
 
Location: Uganda

Join Date: Jan 2015
Posts: 71
Default

if i want to convert lib4seq.sorted.mapped.bam to a fastq file (creating 2 files for paired end) do i need to sort this bam file again?

Last edited by kaps; 05-19-2015 at 02:39 AM. Reason: error correction
kaps is offline   Reply With Quote
Old 05-19-2015, 04:51 AM   #16
kaps
Member
 
Location: Uganda

Join Date: Jan 2015
Posts: 71
Default

I have generated fastq files (R1& R2) from the mapped bam file and I am now interested in assembling those reads using velvet. I need guidance to how to run velveth, velvetg and viewing the output.
kaps is offline   Reply With Quote
Old 05-22-2015, 12:34 AM   #17
kaps
Member
 
Location: Uganda

Join Date: Jan 2015
Posts: 71
Default

Quote:
Originally Posted by Michael.Ante View Post
Hi Kaps,

you should have a look at the Bowtie2 homepage. There, it is explained in detail how the programs work. At the end of the manual is a "Lambda phage example", which has quite an overlap to your problem. It also has a SAMtools downstream section...

Cheers,
Michael
I have used a reference (comprising of 2 sequences varying in length) to do alignment on paired end fastq files. I have generated the sam file but when i try converting it (sam) to bam, it is not working!
here is the feedback: samtools view -bS 4_S4_paired.sam 4_S4_paired.bam
[main_samview] random alignment retrieval only works for indexed BAM or CRAM files.
▒BC▒m▒▒N▒0D▒ܜK▒X=p[▒qU|
▒ ▒▒▒5r[SY▒I▒▒▒O▒sH[@▒▒▒f▒▒▒a0▒w3▒Rʄq▒▒G▒[▒2▒▒Q▒Lk▒֋▒ vU▒1)J▒ei▒M&"▒qq5▒Q^ݒ▒▒▒▒▒`▒ U▒▒▒Q▒ vI▒▒▒f▒▒q▒▒:;▒=▒▒P▒▒U▒▒"n\o3▒ε▒
▒;*(n▒J)*▒▒▒▒֙%▒zM19H▒I▒▒▒▒▒5{▒>▒S▒H▒Qt▒~▒▒▒
▒▒y▒▒

when with -bT, I get a similar feedback as;
samtools view -bT cp4genome.fa 4_S4_paired.sam 4_S4_paired.bam
[samfaipath] build FASTA index...
[fai_build_core] different line length in sequence '4_S4_contig_23'.
[samfaipath] fail to build FASTA index.
[main_samview] random alignment retrieval only works for indexed BAM or CRAM files.
▒BC▒m▒▒N▒0D▒ܜK▒X=p[▒qU|
▒ ▒▒▒5r[SY▒I▒▒▒O▒sH[@▒▒▒f▒▒▒a0▒w3▒Rʄq▒▒G▒[▒2▒▒Q▒Lk▒֋▒ vU▒1)J▒ei▒M&"▒qq5▒Q^ݒ▒▒▒▒▒`▒ U▒▒▒Q▒ vI▒▒▒f▒▒q▒▒:;▒=▒▒P▒▒U▒▒"n\o3▒ε▒
▒;*(n▒J)*▒▒▒▒֙%▒zM19H▒I▒▒▒▒▒5{▒>▒S▒H▒Qt▒~▒▒▒
where is the problem?
kaps is offline   Reply With Quote
Old 05-22-2015, 02:37 AM   #18
Michael.Ante
Senior Member
 
Location: Vienna

Join Date: Oct 2011
Posts: 121
Default

Hi Kaps,

try
Code:
samtools view -S 4_S4_paired.sam | head
if you don't get any problems there, proceed with
Code:
 samtools view -bSh 4_S4_paired.sam > 4_S4_paired.bam
The missing '>' was then the problem. If the sam file doesn't have any header lines, you'll need to provide these additionally.

Cheers,

Michael
Michael.Ante is offline   Reply With Quote
Old 05-22-2015, 05:07 AM   #19
kaps
Member
 
Location: Uganda

Join Date: Jan 2015
Posts: 71
Default

its working.

Thanks
kaps 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 08:37 AM.


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