Question: Getting coding sequence length
0
gravatar for A
4 weeks ago by
A3.8k
A3.8k wrote:

Hi

I have a list of genes

Where I can get the coding sequence length of each gene (in nucleotide) ?

I have a file having gene length but I guess that IS FOR whole gene length not coding sequence part of a gene

Can you help please?

cds gene ncbi • 159 views
ADD COMMENTlink modified 4 weeks ago • written 4 weeks ago by A3.8k
1

You can also use ensemble biomart directly. For Grch37 there is an archive ensembl biomart site http://grch37.ensembl.org/biomart/martview .

  • Select Human genes
  • In attributes -> Structures -> gene -> CDS Length Radio button

The advantage being other metadata can also be pulled in like gene name etc.

ADD REPLYlink written 4 weeks ago by microfuge1.8k

without the actual CDS coordinates ( structural annotations ) , it will be quite hard .

Don't you have a gff or such file with the annotations for all CDSs?

ADD REPLYlink modified 4 weeks ago • written 4 weeks ago by lieven.sterck8.3k

No I don't

You meant I must extract that from gif file?

ADD REPLYlink written 4 weeks ago by A3.8k

if the species is annotated, you could get the CDS in fasta and count the number of bases there for your genes of interest.

ADD REPLYlink written 4 weeks ago by husensofteng260

This code suppose to get CDS length for human but I dob't know why I get error in second line

ens<-read.table("Homo_sapiens.GRCh38.84.gtf/Homo_sapiens.GRCh38.84.gtf",sep="\t",skip=3)
nn1<-which(ens[,3]=="CDS")
genes<-paste0("ENSG",gsub(".*ENSG","",as.character.factor(ens[nn1,9])))
genes<-gsub(";.*","",genes)
transcr<-paste0("ENST",gsub(".*ENST","",as.character.factor(ens[nn1,9])))
transcr<-gsub(";.*","",transcr)
len<-ens[nn1,5]-ens[nn1,4]
df<-cbind.data.frame(genes,transcr,len)
df1<-aggregate(df[,3],by=list(genes,transcr),FUN="sum")
write.csv(df1,"gene,transcript,CDS_length_ list,ens84,grch38.csv")

Getting error here

>genes<paste0("ENSG",gsub(".*ENSG","",as.character.factor(ens[nn1,9])))

Error in as.character.factor(ens[nn1, 9]) :   attempting to coerce non-factor
ADD REPLYlink written 4 weeks ago by A3.8k

First off, you don't need to skip the comment lines in the GTF file, read.table automatically skips lines that begin with comment.char, whose default value is #.

Also, with GTF file, you're better off using data.table and fread than data.frame and read.table.

Unless you're sure ens[nn1, 9] is a factor, don't use as.character.factor. Simply use as.character and let R figure out which method to dispatch to.

ADD REPLYlink written 4 weeks ago by RamRS28k

Have you looked at: A: How to obtain the length of coding regions for the list of genes?

ADD REPLYlink written 4 weeks ago by genomax87k

My genome is GRCH37, does this differ in CDS length?

ADD REPLYlink written 4 weeks ago by A3.8k
1

Get GTF file for GRCh37 here.

zcat gencode.v34lift37.annotation.gtf.gz | awk '{if($3=="exon") print $10"\t"$5-$4}' | sed -e 's/"//g' -e 's/;//' | bedtools groupby -i - -g 1 -c 2 -o sum > Exon_lengths.txt
ADD REPLYlink written 4 weeks ago by genomax87k
1

A, please stop adding answers. Add comments/comment-replies instead.

ADD REPLYlink written 4 weeks ago by RamRS28k
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: 860 users visited in the last hour