SEQanswers

Go Back   SEQanswers > Applications Forums > RNA Sequencing



Similar Threads
Thread Thread Starter Forum Replies Last Post
GAGE Pathview and gene counts JustinH RNA Sequencing 1 07-24-2015 03:28 PM
samilar gene ID in counts file shirley47162928 Bioinformatics 5 06-29-2015 05:00 PM
HTSeq exon counts to gene level Minty RNA Sequencing 10 06-18-2014 04:48 AM
edgeR mean normalized reads counts for each gene antoza Introductions 2 02-11-2014 08:22 AM
DEXSeq gene level counts Julien Roux Bioinformatics 3 11-28-2012 12:31 AM

Reply
 
Thread Tools
Old 04-27-2018, 07:33 AM   #1
alekzs
Junior Member
 
Location: CT

Join Date: Jan 2018
Posts: 7
Default ERCC - no gene counts

Hi!
I have a problem with recovering added ERCC to my RNAseq samples. Briefly, I'm doing smart-seq2 on single human T cells and align with STAR to hg19+ERCC sequences. Both FASTA and GTF file have the ERCCs but when I look at gene count results they don't show up at all, regardless if I use STARs genecount option, HTseq-count or RSEM.
Let's assume I "forgot" to add the ERCC spike-ins and the count is actually 0... Shouldn't the gene names still appear in the downstream files but the value just be 0? If I run Samtools idxstats on the STAR output (sorted or unsorted bam file), it shows the ERCC "chromosomes".
I'm confused by all of this and can't even figure out where in the pipeline my mistake might be. Help!

Here are my commands:
STAR --runMode genomeGenerate --runThreadN 8 --genomeDir indices/STAR --genomeFastaFiles path/to/genomeE.fa --sjdbGTFfile path/to/genesE.gtf --genomeChrBinNbits 12

STAR --runMode alignReads \
--genomeLoad NoSharedMemory \
--genomeDir indices/STAR \
--readFilesIn XX_R1_001.fastq.gz XX_R2_001.fastq.gz \
--outFileNamePrefix /results/ercc/$i \
--quantMode GeneCounts \
--twopassMode Basic \
--outSAMtype BAM Unsorted SortedByCoordinate \
--readFilesCommand zcat

htseq-count --mode=union --idattr=gene_name -f bam -order pos --stranded=no XX-Aligned.bam /path/to/genesE.gtf > XX-gene.count
### I tried stranded=yes or reverse but that didn't help either.

Any pointers highly appreciated!!

Alex
alekzs is offline   Reply With Quote
Old 04-27-2018, 07:48 AM   #2
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,759
Default

You made a "new" reference by appending the fasta ERCC sequences to end of human genome and then created the STAR indexes from this hybrid file?
GenoMax is offline   Reply With Quote
Old 04-27-2018, 09:07 AM   #3
alekzs
Junior Member
 
Location: CT

Join Date: Jan 2018
Posts: 7
Default

Quote:
Originally Posted by GenoMax View Post
You made a "new" reference by appending the fasta ERCC sequences to end of human genome and then created the STAR indexes from this hybrid file?
Yes, I added both FASTA and GTF annotations und used the hybrid!
alekzs is offline   Reply With Quote
Old 04-27-2018, 09:12 AM   #4
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,759
Default

Then I am inclined to speculate that someone forgot to spike the ERCC aliquots. Unless alignments are not being reported since they fail STAR's multi-mapping threshold. Look into that as well.

Did you make the libraries (and add ERCC)?
GenoMax is offline   Reply With Quote
Old 04-27-2018, 09:25 AM   #5
alekzs
Junior Member
 
Location: CT

Join Date: Jan 2018
Posts: 7
Default

Quote:
Originally Posted by GenoMax View Post
Then I am inclined to speculate that someone forgot to spike the ERCC aliquots. Unless alignments are not being reported since they fail STAR's multi-mapping threshold. Look into that as well.

Did you make the libraries (and add ERCC)?
I did everything myself so chances are 50-50 I guess.
Anyhow, even if I didn't add the spike ins, shouldn't the gene names from the reference appear in a gene count file? Like, normal genes get 0 alignments/counts but they're still in the list, right?
alekzs is offline   Reply With Quote
Old 04-27-2018, 10:44 AM   #6
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,759
Default

When you added them to the GTF file they were in the correct format?

Are you able to see alignments for them in the BAM file?
GenoMax is offline   Reply With Quote
Old 04-27-2018, 11:30 AM   #7
alekzs
Junior Member
 
Location: CT

Join Date: Jan 2018
Posts: 7
Default

Quote:
Originally Posted by GenoMax View Post
When you added them to the GTF file they were in the correct format?

Are you able to see alignments for them in the BAM file?
Code:
Tail of FASTA file:
>ERCC-00171 DQ854994 Ac03459967_a1 Ac03460063_a1
CTGGAGATTGTCTCGTACGGTTAAGAGCCTCCGCCCGTCTCTGGGACTATGGACGGGCACGCTCATATCAGGCTATATTTGGTCCGGGTTATTATCGTCGCGGTTACCGTAATACTTCAGATCAGTTAAGTAGGGCCATATGCCTCGGGAATAAGCTGACGGTGACAAGGTTTCCCCCTAATCGAGACGCTGCAATAACACAGGGGCATACAGTAACCAGGCAAGAGTTCAATCGCTTAGTTTCGTGGCGGGATTTGAGGAAAACTGCGACTGTTCTTTAACCAAACATCCGTGCGATTCGTGCCACTCGTAGACGGCATCTCACAGTCACTGAAGGCTATTAAAGAGTTAGCACCCACCATTGGATGAAGCCCAGGATAAGTGACCCCCCCGGACCTTGGAGTTTCATGCTAATCAAAGAAGAGCTAATCCGACGTAAAGTTGCGGCGTTGATTACGCAGGATTGCGACCAAAGAACGAGAAAAAAAAAAAAAAAAAAAAAAAA

Tail of GTF file
>ERCC-00171	ercc	gene	1	506	.	+	.	gene_id "GERCC-00171"; gene_version "1"; gene_name "ERCC-00171"; gene_source "ercc"; gene_biotype "ercc";

samtools view -h 10BTreg02_S290_L003Aligned.sortedByCoord.out.bam ERCC-00171
>NS500597:113:HH5HKBGX5:3:11406:4418:20117	83	ERCC-00171	441	255	9S29M	=	60	-410	ACGACGTAGGTTGCGGCGTTGATTACGCAGGATTGCGA	EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAAAA	NH:i:1	HI:i:1	AS:i:65	nM:i:0
NS500597:113:HH5HKBGX5:3:21612:3414:12240	89	ERCC-00171	442	255	36M	*	0	0	TTGCGGCGTTGATTACGCAGGATTGCTACCAAAGAA	EAEEEEEEAEAEEEEE/EEEEEEEAEEEEEEAAAAA	NH:i:1	HI:i:1	AS:i:33	nM:i:1
(there's many more lines, finds other ERCC numbers as well)

tail -n 20 10BTreg02_S290_L003ReadsPerGene.out.tab
ENSG00000224240	0	0	0
ENSG00000227629	0	0	0
ENSG00000237917	0	0	0
ENSG00000231514	0	0	0
ENSG00000235857	0	0	0
That's all I have to offer.

Last edited by GenoMax; 04-27-2018 at 11:52 AM. Reason: Added [code] tags
alekzs is offline   Reply With Quote
Old 04-27-2018, 11:31 AM   #8
r.rosati
Member
 
Location: Brazil

Join Date: Aug 2015
Posts: 67
Default

Here I am with makeshift solutions, but if you make the BAM into a SAM, you can `grep` it to see if the sequences are there.
r.rosati is offline   Reply With Quote
Old 04-27-2018, 11:43 AM   #9
alekzs
Junior Member
 
Location: CT

Join Date: Jan 2018
Posts: 7
Default

Quote:
Originally Posted by r.rosati View Post
Here I am with makeshift solutions, but if you make the BAM into a SAM, you can `grep` it to see if the sequences are there.
ha, that approach was far easier...

grep "ERCC-" 10B02_3.sam -c
6770

So, yes... they are there, just don't end up in any count file.
alekzs is offline   Reply With Quote
Old 04-27-2018, 11:53 AM   #10
r.rosati
Member
 
Location: Brazil

Join Date: Aug 2015
Posts: 67
Default

I meant like grepping for
CTGGAGATTGTCTCGTACGGTTAAGAGCCTCCGCCC
(or any other fragment in the ERCC controls, I copy-pasted the one you wrote in a previous post)
r.rosati is offline   Reply With Quote
Old 04-27-2018, 11:58 AM   #11
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,759
Default

Can you try featureCounts to do the counts? It will not count multi-mapping reads by default.
GenoMax is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 11:58 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2018, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO