Question: Extracting reads mapping to CDS regions but not the first 45bp after ATG
0
gravatar for Adrian Pelin
7 days ago by
Adrian Pelin2.3k
Canada
Adrian Pelin2.3k wrote:

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 bam bed • 92 views
ADD COMMENTlink modified 7 days ago by Pierre Lindenbaum122k • written 7 days ago by Adrian Pelin2.3k
2
gravatar for Pierre Lindenbaum
7 days ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum122k wrote:

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
(...)
ADD COMMENTlink written 7 days ago by Pierre Lindenbaum122k

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.

ADD REPLYlink written 7 days ago by Adrian Pelin2.3k
1

ENST00000442639 is discarded because there is no stop codon (Incomplete CDS) http://www.ensembl.org/Homo_sapiens/Transcript/Summary?g=ENSG00000134321;r=2:6865828-6887053;t=ENST00000442639

ADD REPLYlink written 7 days ago by Pierre Lindenbaum122k

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.

ADD REPLYlink written 7 days ago by Adrian Pelin2.3k
1

also remove the last 5 codons

should be something like

positions.subList(45,positions.size()-15);

for + strand. The function sublist is https://docs.oracle.com/javase/8/docs/api/java/util/List.html#subList-int-int-

ADD REPLYlink written 7 days ago by Pierre Lindenbaum122k

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

ADD REPLYlink written 7 days ago by Adrian Pelin2.3k

yep ! :-)

ADD REPLYlink written 6 days ago by Pierre Lindenbaum122k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1198 users visited in the last hour