This is a test version of Biostars. For the public version, visit https://www.biostars.org.
no output after running samtools depth

I want to ask another question related to my previous post here: C: compute depth of coverage for certain regions of the genome using samtools depth

I have downloaded bam files of 1000 genomes project phase 3, "not the entire genome but only for certain regions" by using samtools view. I want to compute depth of coverage for these specific regions which can be done by using samtools depth -r .

According to my previous post, samtools depth -r will not work and first I need to reindex my bam files. I re-index them by using 

wget -q -O - "ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/current.tree" | cut -f 1 | grep '.bam$' | grep -Ev -i ".chrom20.|.chrom11." | parallel --dry-run /data/programs/samtools-1.3.1/samtools view -b -o {/} {} "16:32297737-33456560"

however, when I try samtools depth -r after re-indexing, 

/data/programs/samtools-1.3.1/samtools depth -r 16:32297737-33456560 NA21137.mapped.ILLUMINA.bwa.GIH.low_coverage.20120522.bam  > bb

my output file is empty, I also get an error message that

Warning: The index file is older than the data file: NA21137.mapped.ILLUMINA.bwa.GIH.low_coverage.20120522.bam.bai

I wonder what is wrong here? dos anyone have any idea? thanks

samtools bam
• 3,355 views • link
written by Ana Ana

1 answer

Index the bam file and try again. 

wget -q -O - "ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/current.tree" | cut -f 1 | grep '.bam$' | grep -Ev -i ".chrom20.|.chrom11." | parallel --dry-run /data/programs/samtools-1.3.1/samtools view -b -o {/} {} "16:32297737-33456560"

Above script will not fetch any bam as you are dry running parallel. To download you are supposed to remove --dry-run.

• 0 views • link
written by cpad0112 cpad0112

thanks @cpad0112, unfortunately I really do not understand your response and what do you mean by indexing bam files and trying the code again! I downloaded the bam files by this command

wget -q -O - "ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/current.tree" | cut -f 1 | grep '.bam$' | while read B; do echo -n "$B " && /data/programs/samtools-1.3.1/samtools  view -b -o ${B##*/} "http://ftp.1000genomes.ebi.ac.uk/vol1/$B" "16:32297737-33456560"  && rm *.chrom20.* && rm *.chrom11.* ; done

the outputs are bam and bam.bai files (bam.bai are index files as far as I know). later on I re-index the files by this

wget -q -O - "ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/current.tree" | cut -f 1 | grep '.bam$' | grep -Ev -i ".chrom20.|.chrom11." | parallel --dry-run /data/programs/samtools-1.3.1/samtools view -b -o {/} {} "16:32297737-33456560"

Can you explain it bit more for me, I am totally new in human genome project, a bit confused by what you mean exactly.

• 0 views • link
written by Ana Ana

wget -q -O - "ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/current.tree" | cut -f 1 | grep '.bam$' | grep -Ev -i ".chrom20.|.chrom11." | parallel --dry-run /data/programs/samtools-1.3.1/samtools view -b -o {/} {} "16:32297737-33456560" is not for indexing the bam files. This a different way to download the bam files from internet. samtools index input.bam is the right command to index bam files.

Following is a test run for indexing bam and calculating depth by samtools on a single file, in a temporary directory. If these steps are successful, then you can apply them to all the bams in the directory.

  1. make a new directory (name it like test)
  2. Copy NA21137.mapped.ILLUMINA.bwa.GIH.low_coverage.20120522.bam and NA21137.mapped.ILLUMINA.bwa.GIH.low_coverage.20120522.bam.bai (copy, don't move the file)
  3. cd to the new directory
  4. You should see two files: with .bam and .bam.bai extension.
  5. Delete .bai file (rm NA21137.mapped.ILLUMINA.bwa.GIH.low_coverage.20120522.bam.bai)
  6. Now you should see only bam file

  7. Index the bam file with following command:

    samtools index NA21137.mapped.ILLUMINA.bwa.GIH.low_coverage.20120522.bam

  8. Now you should see two files with .bam and .bam.bai extensions

  9. Run samtools depth (example given below):

    samtools depth -r 16:32297737-32297837 NA21137.mapped.ILLUMINA.bwa.GIH.low_coverage.20120522.bam

  10. Step 9 would print the result to the screen. Result can be saved in a separate file.

samtools depth -r 16:32297737-32297837 NA21137.mapped.ILLUMINA.bwa.GIH.low_coverage.20120522.bam > NA21137.mapped.ILLUMINA.bwa.GIH.low_coverage.20120522. text.

Note: I edited post mentioning parallel that parallel is used for downloading bams, as alternative to bash script provided in that post.

• 0 views • link
written by cpad0112 cpad0112

If above steps are successful, navigate to wherever bam files are. Try running following command:

$ for i in *.bam;do  samtools index $i;done

This would index all the bam files. Try running samtools depth again on bam files.

• 0 views • link
written by cpad0112 cpad0112

thanks so much, excellent help

• 0 views • link
written by Ana Ana

Log in to answer this question.

More posts like this