1
$\begingroup$

I have a RNA aligned file without any header info. I need to calculate the coverage of each mapped base for per chromosome. Is it possible to do it with perl ? 

NS500377:28:H16A7BGXX:2:21205:16038:7658        99      chr2   3000482 255     46M     =       3000502 56      GTCCTACAGTGTGCATTTCTCGTTTTCCACGTTTCT    AAAA)FFFFFFFFFAFFFFFFFFFFF.F.FFFFFAF    XA:i:2  MD:Z:21T4T7T1   NM:i:3
NS500377:28:H16A7BGXX:3:11503:24953:10427       99      chr2   3000482 255     50M     =       3000513 71      GTCCTACAGTGTGCATTTCTCATTTTTCAAGTTTTTCAGT        AAAAAFFFFFFFFFFFFAFFFFFF7AAFFFFFFFFFFFFF
        XA:i:1  MD:Z:21T7C10    NM:i:2
NS500377:28:H16A7BGXX:3:12404:23929:9392        163     chr2   3000482 255     49M     =       3000483 40      GTCCTACAGTGTGCATTTCCCTTTTTTCACGTTTTTTAG AAFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF XA:i:1  MD:Z:19T16C2    NM:i:2
NS500377:28:H16A7BGXX:3:12404:23929:9392        83      chr2   3000483 255     49M     =       3000482 -40     TCCTACAGTGTGCATTTCCCTTTTTTCACGTTTTTTAGT AFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF XA:i:2  MD:Z:18T16C3    NM:i:2
NS500377:28:H16A7BGXX:3:21606:15675:1639        83      chr2   3000483 255     50M     =       3000480 -43     TCTTACAGTGTGCATTTCTCTTTTTTCACGTTTTTCAGTG        AFFFFFFFFFAFFFFFFFFFFFFFFFFFFFFFFFFAAAAA
        XA:i:0  MD:Z:2C37       NM:i:1
NS500377:28:H16A7BGXX:4:12605:18505:11088       83      chr2   3000486 255     50M     =       3000481 -45     TACAATGTGCATTTCTCTTTTTTCACGTTTTTCAGTATTT        FFFFFFFAFFF)FFFFFFFFF<FFFF<AFFFFFFFAA<AA
        XA:i:1  MD:Z:4G31G3     NM:i:2
NS500377:28:H16A7BGXX:2:13303:20634:12085       99      chr2   3000487 255     49M     =       3000521 73      ACAGTGTACATTTCTCATTTTTCACGTTTTTCAGTGTTA AAAAAFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF. XA:i:2  MD:Z:7G8T21T0   NM:i:3
Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

You are looking for bedtools genomecov, specifically -d option. See https://bedtools.readthedocs.io/en/latest/content/tools/genomecov.html , for example:

bedtools genomecov -i A.bed -g my.genome -d

I recommend not reinventing the wheel: use bedtools instead of reimplementing the same functionality in Perl.

Improve this answer
$\endgroup$
3
  • $\begingroup$ thanks, but bedtools doesn't work in this case as there is no header info. $\endgroup$
    – Gene_of_Thrones
    Jan 28, 2020 at 14:21
  • $\begingroup$ Can you add it from another bam file using samtools? Or add a fake header? $\endgroup$
    – Timur Shtatland
    Jan 28, 2020 at 14:25
  • $\begingroup$ i did add a header from a different sample to calculate the coverage... just wondering if we can do it using awk. $\endgroup$
    – Gene_of_Thrones
    Feb 10, 2020 at 13:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.