This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Calculating amount of reads on certain chromosomes in bam file

Hi,

I'm looking for a nice way to calculate the percentage of reads on specifique chromosomes (MT und (1-23,X,Y) und unplaced scaffold).

I know I can get the read counts per chromosome with samtools idxstats, however I have a lot of bam files and I would like to automate the calculation. The problem is, I'm struggling with basic batch text maniplulation and would therefore appreciate any help you can give me (specific or general).

Edit: I forgot to mention, that I have multiple bam files, which are indexed and only consist of uniquely mapped reads.

rna-seq bam samtools
• 6,204 views • link
updated by WouterDeCoster WouterDeCoster • written by caggtaagtat caggtaagtat

Do you have some familiarity with R or python?

• 0 views • link
written by Sean Davis Sean Davis

Yes, I always work with R

• 0 views • link
written by caggtaagtat caggtaagtat

So, can you update your question with where, specifically, you are getting stuck with processing many files with samtools idxstats? Ideally, post any code you have tried.

• 0 views • link
written by Sean Davis Sean Davis

I will try to do it with R now. I just thought, there might be a short command in bash :)

• 0 views • link
written by caggtaagtat caggtaagtat

You can certainly loop over your files in bash with something like (untested):

for i in `ls *bai`; do j=`echo $i | sed s/\.bai//`; samtools idxstats $i > $j.idxstats; done

I suspect you will still want to "analyze" the data, so R will likely be involved at some point.

• 0 views • link
written by Sean Davis Sean Davis

2 answers

General workflow you could try given that you use R.

  1. Look at using dir() to get a list of index files.
  2. Look at using sprintf() to create a command string to run samtools idxstats given a filename (which you have from #1), piping the results to a unique file for each sample.
  3. Look at using system() to run the command string from 2.
  4. Look at using a for loop loop to run #3 over your list of index files from #1.
  5. Define a function to read a single idxstats output file.
  6. Loop using sapply or lapply over all the idxstat output files
  7. Combine output of 6 into a single matrix, with samples in columns and chromosomes in rows.
• 0 views • link
written by Sean Davis Sean Davis

Thank you, I will do that

• 0 views • link
written by caggtaagtat caggtaagtat

I'd suggest to also having a look at mosdepth and indexcov.

• 0 views • link
written by WouterDeCoster WouterDeCoster

You have the same link for both programs mentioned.

• 0 views • link
written by GenoMax GenoMax

Log in to answer this question.

More posts like this