Hi guys, I'm a newbie in the RNA Seq world. I'm analyzing some RNA seq data (HiSeq) and doing this, I'm following the best practice guideline available online form different sources. Briefly I performed the alignment to the reference genome, the merge, I removed duplicates, then I performed the sorting, indexing and finally using bedtools multicov I ended up having the reads count per RefSeq. My point is that I would like to have a gene level expression measure but after the counts I have reads counts per transcript variants. in other words I have multiple transcripts per gene with corresponding reads count but I would like to have one reads counts per gene. How do you deal with this issue?
Thanks in advance
B.
•
written
by
Bnf83
1 answer
Firstly I would agree that removing duplicates is generally only something that should done if there are reasons, and not as a default.
It sounds like you used a reference set of refseq transcripts when counting with multicov. As far as I can see multicov performs a naive count of the reads in each interval given, which means the counts you observe on transcript level are not really meaningful as transcripts have many overlaps.
To get what you want you simply need to count features on a gene level. Both featureCounts as mentioned before and for example htseq-count can be used for this, as they can count reads on a gene metalevel which will count reads from all possible exons in each gene and summarize the results. I am unsure about the annotation quality of other species, but in human I have usually used the lastest version of gencode annotation in gtf or gff format (which includes gene and exon information) together with htseq-count/featureCounts. Note that your annotation version and reference genome used in alignment must match. You just have to make sure the meta feature you are counting is gene, and the feature is exons.
Having looked recently at some human RiboSeq data, I'm not sure how accurate the biotypes tend to be. I've seen a good number of non-protein coding getting translated.
The organism won't matter.
Log in to answer this question.
-
read count for RNA sequencing via nanopore
written by gat •
-
After read mapping, count exons or genes?
written by tianshenbio
-
Millions of Reads Needed for Differential Gene Expression
written by contact •
-
Chip seq (ac) reads count per pick
written by Bnf83
-
CDS read counts per gene
written by ashkan
-
Do you use anti-sens reads for differential gene expression analysis with Stranded Paired-End Libra…
written by ZheFrench
-
RNA Seq coverage
written by Bnf83
-
How to define the gene length for RPKM calculation
written by Bnf83
-
bedtools multicov not running
written by Bnf83
-
Are indel realignment and base recalibration necessary for RNAseq data?
written by Bnf83
What do you want to use the counts for? Normally you'd use featureCounts for that (and not remove duplicates, which is generally a bad idea in RNAseq).
I only would like to know if the genes are differentially expressed. I don't want to know if a specific splicing variant is differentially expressed. Only the gene.
As Devon wrote, use featureCounts (or htseq-count or ...), which you can use to count per gene.