This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Getting coding sequence length

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?

ncbi gene cds
• 5,839 views • link
updated by newbio17 newbio17 • written by zizigolu zizigolu

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?

• 0 views • link
written by lieven.sterck lieven.sterck

No I don't

You meant I must extract that from gif file?

• 0 views • link
written by zizigolu zizigolu

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

• 0 views • link
written by husensofteng husensofteng

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
• 0 views • link
written by zizigolu zizigolu

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.

• 0 views • link
written by Ram Ram

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

• 0 views • link
written by zizigolu zizigolu

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
• 0 views • link
written by GenoMax GenoMax

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.

• 0 views • link
written by newbio17 newbio17

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

• 0 views • link
written by Ram Ram

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.

• 0 views • link
written by microfuge microfuge

0 answers

No answers yet.

Log in to answer this question.

More posts like this