Which feature file for htseq-count for non coding elements of ribodepleted samples?

Jane M · 05-10-2016, 01:52 AM

Dear all,

Sorry if this question has been asked before, I have not found similar topics on the web, but there are probably...
Until now, I used a gtf downloaded from UCSC table (genome: human, assembly: hg19, group:genes and gene predictions, track: RefSeq genes, table: refFlat) as a feature file, with gene name, in GFF format for counting reads in genes with htseq-count.
I have one experiment with ribodepleted samples (TruSeq total RNA Stranded). I first counted the reads and performed differential expression analysis using such a gtf. In this gtf downloaded few weeks ago, there are 26688 genes, including 913 "LINC*". This number of lincRNA seems low, so I wonder if this table is comprehensive for ribodepleted experiment.
Can you please tell me which reference you use when interesting in non coding elements of ribodepleted experiments?

Thank you for your feedback,
Jane

dpryan · 05-10-2016, 11:58 AM

UCSC references...kind of suck. Use Gencode/Ensembl and you'll get more complete coding and non-coding transcripts. Having said that, for lincRNAs you might want to check out lincrnadb or RNAcentral.

hyeonjipark · 05-10-2016, 10:21 PM

I used gencode. but no counts.
perhaps what column used htseq-count ?

dpryan · 05-11-2016, 12:18 AM

Presumably the lack of counts is due to the difference in chromosome names.

Jane M · 05-11-2016, 01:21 AM

Thank you dpryan for your answer.

Quote:

Originally Posted by dpryan

UCSC references...kind of suck.

Don't you thing this refFlat reference is sufficient for experiments with polyA RNA selection?

Quote:

Originally Posted by dpryan

Use Gencode/Ensembl and you'll get more complete coding and non-coding transcripts. Having said that, for lincRNAs you might want to check out lincrnadb or RNAcentral.

Since I am working with hg19 annotation, I downloaded the gtf from http://www.gencodegenes.org/releases/19.html for "ALL" regions, that is the gencode.v19.chr_patch_hapl_scaff.annotation file. I am not familiar with these annotation files. Is it the right one?

I am currently looking at what this file contains: number of genes (name, ID, status), transcripts, ...
It has not the same format like the refFlat file, so I expect some issues for read counting.

dpryan · 05-11-2016, 01:23 AM

The UCSC file annotations are never that good. The file you downloaded is fine, though again the chromosome names probably differ, which is what's causing problems. I think featureCounts can handle the chromosome naming difference internally, so try using that instead (it's much faster anyway).

Jane M · 05-11-2016, 02:46 AM

Quote:

Originally Posted by dpryan

The UCSC file annotations are never that good. The file you downloaded is fine, though again the chromosome names probably differ, which is what's causing problems. I think featureCounts can handle the chromosome naming difference internally, so try using that instead (it's much faster anyway).

Thank you again. No problem yet, since I did not try to count. Looking at the file: now, I have ~56600 different gene names, compared to ~26600 with refFlat. The file contains twice more lines. Looking forward to see how the number of LINC and AS will increase.
chr1->22, chrX, chrY have the same names, but the additional chromosomes have indeed different names: chr17_ctg5_hap1, chr17_gl000205_random in refFlat and GL000191.1, GL000192.1. I start the tests right now!

Jane M · 05-11-2016, 07:17 AM

Compared to how I ran htseq-count previously, I added -i option to get counts for each gene_name and not gene_ID. I could run htseq count on this new gtf without any problem

Code:

htseq-count --format=bam --order=name --mode=union -a 20 -i gene_name --stranded=reverse mybam.bam gencode.v19.chr_patch_hapl_scaff.annotation.gtf > mynewcounts.gencod.htseqCount

But I still have a couple of questions:
1. With the genecode annotation, we can find how many genes are attributed to each gene_type (30 gene types), like protein coding, lincRNA, pseudogene, rRNA,... Do you know where to find this information for the refFlat annotation?

2. I am a bit confused by gene_name and gene_ID:
When counting in gencode.v19.chr_patch_hapl_scaff.annotation.gtf the number of genes ("gene" in field 3), I get 63,568:

Code:

cat gencode.v19.chr_patch_hapl_scaff.annotation.gtf | cut -f3  | wc -l
63568

but the number of unique gene_name is 56,629:

Code:

cat gencode.v19.chr_patch_hapl_scaff.annotation.gtf | cut -d';' -f5 | sort | uniq | wc -l
56629

In the output file of htseq count, 56,629 are listed.

Does someone know what are those genes with same gene_name but different gene_ID?

Code:

cat gencode.v19.chr_patch_hapl_scaff.annotation.gene.gtf | cut -f9  | sort | uniq | wc -l
63568

3. Finally, I noticed that the counting changed a lot between the 2 annotations, even for well known genes as TP53, the number of reads attributed to this gene - and others - double! I am surprised that the annotation of such very well known genes changes. Is it normal?

Any clarification would be greatly appreciated.

dpryan · 05-11-2016, 01:04 PM

1. Use things like "cut" and "uniq" to determine this. This isn't something you need to look up, just determine it yourself.
2. How does one define a gene? Is it a location, a sequence, something else? If you have essentially the same sequence on different chromosomes and both are expressed are they the same gene or different ones? In such cases, gencode/ensembl will give each instance a unique ID. UCSC will give each instance the same ID in such cases, which is a good way to completely break a LOT of programs.This is why one should normally quantify by gene ID. You can add gene names after everything is analysed.
3. UCSC annotations are rather minimalistic.

Jane M · 05-12-2016, 06:35 AM

Quote:

Originally Posted by dpryan

1. Use things like "cut" and "uniq" to determine this. This isn't something you need to look up, just determine it yourself.

Well, there is no description in the tables I downloaded from UCSC. Otherwise, I could indeed check as I did with Genecode annotation. Here are some file header:

refFlat file (RefSeq), that I used until now:

Code:

chr1	hg19_refFlat	exon	11874	12227	0.000000	+	.	gene_id "DDX11L1"; transcript_id "DDX11L1"; 
chr1	hg19_refFlat	exon	12613	12721	0.000000	+	.	gene_id "DDX11L1"; transcript_id "DDX11L1"; 
chr1	hg19_refFlat	exon	13221	14409	0.000000	+	.	gene_id "DDX11L1"; transcript_id "DDX11L1"; 
chr1	hg19_refFlat	exon	14362	14829	0.000000	-	.	gene_id "WASH7P"; transcript_id "WASH7P"; 
chr1	hg19_refFlat	exon	14970	15038	0.000000	-	.	gene_id "WASH7P"; transcript_id "WASH7P"; 
chr1	hg19_refFlat	exon	15796	15947	0.000000	-	.	gene_id "WASH7P"; transcript_id "WASH7P";

refGene file (RefSeq):

Code:

chr1	hg19_refGene	start_codon	67000042	67000044	0.000000	+	.	gene_id "NM_032291"; transcript_id "NM_032291"; 
chr1	hg19_refGene	CDS	67000042	67000051	0.000000	+	0	gene_id "NM_032291"; transcript_id "NM_032291"; 
chr1	hg19_refGene	exon	66999639	67000051	0.000000	+	.	gene_id "NM_032291"; transcript_id "NM_032291"; 
chr1	hg19_refGene	CDS	67091530	67091593	0.000000	+	2	gene_id "NM_032291"; transcript_id "NM_032291"; 
chr1	hg19_refGene	exon	67091530	67091593	0.000000	+	.	gene_id "NM_032291"; transcript_id "NM_032291"; 
chr1	hg19_refGene	CDS	67098753	67098777	0.000000	+	1	gene_id "NM_032291"; transcript_id "NM_032291"; 
chr1	hg19_refGene	exon	67098753	67098777	0.000000	+	.	gene_id "NM_032291"; transcript_id "NM_032291"; 
chr1	hg19_refGene	CDS	67101627	67101698	0.000000	+	0	gene_id "NM_032291"; transcript_id "NM_032291";

knownGenes file (UCSC):

Code:

chr1	hg19_knownGene	exon	11874	12227	0.000000	+	.	gene_id "uc010nxr.1"; transcript_id "uc010nxr.1"; 
chr1	hg19_knownGene	exon	12646	12697	0.000000	+	.	gene_id "uc010nxr.1"; transcript_id "uc010nxr.1"; 
chr1	hg19_knownGene	exon	13221	14409	0.000000	+	.	gene_id "uc010nxr.1"; transcript_id "uc010nxr.1"; 
chr1	hg19_knownGene	start_codon	12190	12192	0.000000	+	.	gene_id "uc010nxq.1"; transcript_id "uc010nxq.1"; 
chr1	hg19_knownGene	CDS	12190	12227	0.000000	+	0	gene_id "uc010nxq.1"; transcript_id "uc010nxq.1"; 
chr1	hg19_knownGene	exon	11874	12227	0.000000	+	.	gene_id "uc010nxq.1"; transcript_id "uc010nxq.1"; 
chr1	hg19_knownGene	CDS	12595	12721	0.000000	+	1	gene_id "uc010nxq.1"; transcript_id "uc010nxq.1"; 
chr1	hg19_knownGene	exon	12595	12721	0.000000	+	.	gene_id "uc010nxq.1"; transcript_id "uc010nxq.1"; 
chr1	hg19_knownGene	CDS	13403	13636	0.000000	+	0	gene_id "uc010nxq.1"; transcript_id "uc010nxq.1"; 
chr1	hg19_knownGene	stop_codon	13637	13639	0.000000	+	.	gene_id "uc010nxq.1"; transcript_id "uc010nxq.1"; 
chr1	hg19_knownGene	exon	13403	14409	0.000000	+	.	gene_id "uc010nxq.1"; transcript_id "uc010nxq.1";

There might exist one file with comprehensive description.

Quote:

Originally Posted by dpryan

2. How does one define a gene? Is it a location, a sequence, something else? If you have essentially the same sequence on different chromosomes and both are expressed are they the same gene or different ones? In such cases, gencode/ensembl will give each instance a unique ID. UCSC will give each instance the same ID in such cases, which is a good way to completely break a LOT of programs.This is why one should normally quantify by gene ID. You can add gene names after everything is analysed.

Thank you for the explanation.

Quote:

Originally Posted by dpryan

3. UCSC annotations are rather minimalistic.

Ok, but I am very surprised for annotation of well known genes.

GenoMax · 05-12-2016, 06:43 AM

Why are you not using the GTF file from UCSC instead of all those other files. You could create that from Table Browser. (iGenomes bundle does not have non-coding genes).

Edit: @Vikas Bansal has a solution to get the non-coding elements (you can choose GTF output in table browser) for hg19 here.

Jane M · 05-12-2016, 07:26 AM

Quote:

Originally Posted by GenoMax

Why are you not using the GTF file from UCSC instead of all those other files. You could create that from Table Browser. (iGenomes bundle does not have non-coding genes).

I was actually using only the refFlat gtf from UCSC. It looks like what I showed.
I added the refGene and knownGenes gtf files to show that these files do not contain description neither.
There seem to be some non coding elements in refFlat gtf: ~900 LINC, ...

Quote:

Originally Posted by GenoMax

Edit: @Vikas Bansal has a solution to get the non-coding elements (you can choose GTF output in table browser) for hg19 here.

Thank you, I will take a look. But it is probably easier to use genecode to get everything from a single gtf file.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
DEXSeq Using Counts File From htseq-count	FuzzyCoder	Bioinformatics	20	01-04-2016 12:18 AM
Error with GTF file when using htseq-count	MDonlin	Bioinformatics	13	01-13-2015 09:29 AM
Isolating non coding RNA using bedtools/samtools/htseq-count	tirohia	Bioinformatics	3	09-30-2014 12:13 AM
HTseq-count feature type choice for RNAseq	shangzhong0619	Bioinformatics	3	08-07-2014 01:29 PM
multiBamCov or htseq-count to count read per feature ?	NicoBxl	Bioinformatics	1	07-03-2012 03:05 AM

05-10-2016, 01:52 AM	#1
Jane M Senior Member Location: Paris Join Date: Aug 2011 Posts: 239	Which feature file for htseq-count for non coding elements of ribodepleted samples? Dear all, Sorry if this question has been asked before, I have not found similar topics on the web, but there are probably... Until now, I used a gtf downloaded from UCSC table (genome: human, assembly: hg19, group:genes and gene predictions, track: RefSeq genes, table: refFlat) as a feature file, with gene name, in GFF format for counting reads in genes with htseq-count. I have one experiment with ribodepleted samples (TruSeq total RNA Stranded). I first counted the reads and performed differential expression analysis using such a gtf. In this gtf downloaded few weeks ago, there are 26688 genes, including 913 "LINC*". This number of lincRNA seems low, so I wonder if this table is comprehensive for ribodepleted experiment. Can you please tell me which reference you use when interesting in non coding elements of ribodepleted experiments? Thank you for your feedback, Jane

05-11-2016, 07:17 AM	#8
Jane M Senior Member Location: Paris Join Date: Aug 2011 Posts: 239	Compared to how I ran htseq-count previously, I added -i option to get counts for each gene_name and not gene_ID. I could run htseq count on this new gtf without any problem Code: htseq-count --format=bam --order=name --mode=union -a 20 -i gene_name --stranded=reverse mybam.bam gencode.v19.chr_patch_hapl_scaff.annotation.gtf > mynewcounts.gencod.htseqCount But I still have a couple of questions: 1. With the genecode annotation, we can find how many genes are attributed to each gene_type (30 gene types), like protein coding, lincRNA, pseudogene, rRNA,... Do you know where to find this information for the refFlat annotation? 2. I am a bit confused by gene_name and gene_ID: When counting in gencode.v19.chr_patch_hapl_scaff.annotation.gtf the number of genes ("gene" in field 3), I get 63,568: Code: cat gencode.v19.chr_patch_hapl_scaff.annotation.gtf \| cut -f3 \| wc -l 63568 but the number of unique gene_name is 56,629: Code: cat gencode.v19.chr_patch_hapl_scaff.annotation.gtf \| cut -d';' -f5 \| sort \| uniq \| wc -l 56629 In the output file of htseq count, 56,629 are listed. Does someone know what are those genes with same gene_name but different gene_ID? Code: cat gencode.v19.chr_patch_hapl_scaff.annotation.gene.gtf \| cut -f9 \| sort \| uniq \| wc -l 63568 3. Finally, I noticed that the counting changed a lot between the 2 annotations, even for well known genes as TP53, the number of reads attributed to this gene - and others - double! I am surprised that the annotation of such very well known genes changes. Is it normal? Any clarification would be greatly appreciated. Last edited by Jane M; 05-11-2016 at 07:49 AM.

05-12-2016, 06:43 AM	#11
GenoMax Senior Member Location: East Coast USA Join Date: Feb 2008 Posts: 7,080	Why are you not using the GTF file from UCSC instead of all those other files. You could create that from Table Browser. (iGenomes bundle does not have non-coding genes). Edit: @Vikas Bansal has a solution to get the non-coding elements (you can choose GTF output in table browser) for hg19 here. Last edited by GenoMax; 05-12-2016 at 06:54 AM.

05-10-2016, 11:58 AM	#2
dpryan Devon Ryan Location: Freiburg, Germany Join Date: Jul 2011 Posts: 3,480	UCSC references...kind of suck. Use Gencode/Ensembl and you'll get more complete coding and non-coding transcripts. Having said that, for lincRNAs you might want to check out lincrnadb or RNAcentral.

05-10-2016, 10:21 PM	#3
hyeonjipark Junior Member Location: korea Join Date: Dec 2015 Posts: 2	I used gencode. but no counts. perhaps what column used htseq-count ?

05-11-2016, 12:18 AM	#4
dpryan Devon Ryan Location: Freiburg, Germany Join Date: Jul 2011 Posts: 3,480	Presumably the lack of counts is due to the difference in chromosome names.

05-11-2016, 01:23 AM	#6
dpryan Devon Ryan Location: Freiburg, Germany Join Date: Jul 2011 Posts: 3,480	The UCSC file annotations are never that good. The file you downloaded is fine, though again the chromosome names probably differ, which is what's causing problems. I think featureCounts can handle the chromosome naming difference internally, so try using that instead (it's much faster anyway).

05-11-2016, 01:04 PM	#9
dpryan Devon Ryan Location: Freiburg, Germany Join Date: Jul 2011 Posts: 3,480	1. Use things like "cut" and "uniq" to determine this. This isn't something you need to look up, just determine it yourself. 2. How does one define a gene? Is it a location, a sequence, something else? If you have essentially the same sequence on different chromosomes and both are expressed are they the same gene or different ones? In such cases, gencode/ensembl will give each instance a unique ID. UCSC will give each instance the same ID in such cases, which is a good way to completely break a LOT of programs.This is why one should normally quantify by gene ID. You can add gene names after everything is analysed. 3. UCSC annotations are rather minimalistic.