Calculating percentage of reads aligning to a subject

[kaps]

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]

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]

Thanks Michael,

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

[GenoMax]

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]

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]

Originally Posted by kaps

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]

Hi, Michael

Thanks for pointing this to me!

[kaps]

Originally Posted by Michael.Ante

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]

Have you created the index with

samtools index lib4seq.sorted.bam

?
If yes how does your bam-dile header looks like?

 samtools view -H lib4seq.sorted.bam

[kaps]

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

Thanks

[kaps]

Originally Posted by Michael.Ante

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]

Originally Posted by Michael.Ante

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]

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

[Michael.Ante]

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

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]

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]

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]

Originally Posted by Michael.Ante

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]

Hi Kaps,

try

samtools view -S 4_S4_paired.sam | head

if you don't get any problems there, proceed with

 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]

its working.

Thanks