This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Extracting reads mapping to CDS regions but not the first 45bp after ATG

I have an RNA-seq (RPFl Ribosomal Profiling) experiment where I want to focus my analysis on reads mapping to CDS of the human genome, but not the first 15 codons (i.e. first 45bp after ATG start codon).

For this I am trying to build a bed file with coordinates of all CDS and then remove the first 45bp after ATG. There are a lot of issues doing this manually, as for a lot of genes, the first CDS features in a GTF file is shorter than 45bp.

Is there a tool that can do this easily? Thanks.

rna-seq bed bam
• 2,640 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by Adrian Pelin Adrian Pelin

1 answer

using BioAlcidaeJdk http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html

(not checked)

$ java -jar dist/bioalcidaejdk.jar -F GTF -f bioalcidae.code  Homo_sapiens.GRCh37.87.gtf.gz 

3   57557938    57558025    ENST00000303436
3   57561274    57561400    ENST00000303436
3   57563042    57563114    ENST00000303436
3   57569624    57569734    ENST00000303436
3   57570110    57570191    ENST00000303436
3   57582800    57582822    ENST00000303436
3   57557938    57558025    ENST00000496292
3   57561274    57561400    ENST00000496292
3   57563042    57563114    ENST00000496292
3   57569624    57569734    ENST00000496292
(...)
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Ok got it to work. It seems that it only picks one transcript per gene. For instance RSAD2 has 5 different transcripts with 2 of these transcripts represented by CDS features (ENST00000382040 and ENST00000442639). However command above only pulls out ENST00000382040. Based on what criteria does it decide which transcript to fish out? Thanks.

• 0 views • link
written by Adrian Pelin Adrian Pelin

Understood, makes sense. Pierre, thanks a lot for your help. I expected a simpler code and doubt I would be able to tailor it for an additional function. Is it difficult to also remove the last 5 codons (15 bp before stop codon) simultaneous to removing the first 45bp? Sorry I didn't ask right away, thought I would figure it out myself.

• 0 views • link
written by Adrian Pelin Adrian Pelin

Oh that's simple! And "positions = positions.subList(15,positions.size()-45);" after "else" indicating negative strand?

• 0 views • link
written by Adrian Pelin Adrian Pelin

Log in to answer this question.

More posts like this