Getting coding sequence length

Question

[zizigolu]

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?

[microfuge]

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.

[lieven.sterck]

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?

[zizigolu]

No I don't

You meant I must extract that from gif file?

[husensofteng]

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

[zizigolu]

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

[Ram]

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.

[GenoMax]

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

[zizigolu]

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

[GenoMax]

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

[newbio17]

Just a side note:

It might be safer to use transcript columns instead of gene columns (i.e. print $12 instead of print $10) since depending on the annotation an individual might be using, the sum of exon lengths could dramatically change if there are multiple transcripts for the same gene. It could be worthwhile to do grep ENSEMBL before awk in case your GENCODE annotation has HAVANA entries as well.

[Ram]

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