SEQanswers

SEQanswers (http://seqanswers.com/forums/index.php)
-   General (http://seqanswers.com/forums/forumdisplay.php?f=16)
-   -   Calculating percentage of reads aligning to a subject (http://seqanswers.com/forums/showthread.php?t=56569)

kaps 04-08-2015 01:01 AM

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?

Michael.Ante 04-08-2015 02:37 AM

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).

kaps 04-11-2015 02:38 AM

Thanks Michael,

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

GenoMax 04-11-2015 03:39 AM

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?

kaps 04-11-2015 04:57 AM

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)

Michael.Ante 04-14-2015 12:23 AM

Quote:

Originally Posted by kaps (Post 169912)
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

kaps 04-19-2015 11:57 PM

Hi, Michael

Thanks for pointing this to me!

kaps 04-22-2015 01:54 AM

Quote:

Originally Posted by Michael.Ante (Post 168961)
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?

Michael.Ante 04-22-2015 02:42 AM

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

kaps 04-22-2015 06:05 AM

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

Thanks

kaps 04-23-2015 05:43 AM

Quote:

Originally Posted by Michael.Ante (Post 168961)
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 05-11-2015 11:04 PM

Quote:

Originally Posted by Michael.Ante (Post 168961)
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?

GenoMax 05-12-2015 03:20 AM

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 :o

Michael.Ante 05-12-2015 04:24 AM

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

kaps 05-19-2015 02:38 AM

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?

kaps 05-19-2015 04:51 AM

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 05-22-2015 12:34 AM

Quote:

Originally Posted by Michael.Ante (Post 170044)
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?

Michael.Ante 05-22-2015 02:37 AM

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

kaps 05-22-2015 05:07 AM

its working.

Thanks


All times are GMT -8. The time now is 06:01 AM.

Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.