Zero output for sel.rn=rowSums(cnts) != 0

Parharn · 06-16-2014, 04:46 AM

I get a problem when I want to try step 2 from "RNA-seq data pathway and gene-set analysis workflow" by Weijun Luo. I get zero rows with value !=0.

Code:

> p.cnts=assay(gnCnt)
> cnts=p.cnts
> dim(cnts)
[1] 7017    2
> sel.rn=rowSums(cnts) !=0
> cnts=cnts[sel.rn,]
> dim(cnts)
[1] 0 2

When I check my BAM files with IGV it looks that I have mapped reads all over the genome.
Any advise?

bigmw · 06-17-2014, 11:30 AM

I haven’t seen your data and exactly how you did your reads mapping and counting. But please make sure on two things:
First, you use the same version of reference genome at the TopHat step and in your TranscriptDb object.
Second, you specify the right parameters in the reads counting step. They are different for single- vs paired-end data. The workflow example is for paired-end data. Single end data are treated differently:
…
flag <- scanBamFlag(isNotPrimaryRead=FALSE, isProperPair=NA)
param <- ScanBamParam(flag=flag)
gnCnt <- summarizeOverlaps(exByGn, bamfls, mode="Union", ignore.strand=TRUE, single.end=TRUE, param=param)
…

Please check help info for details:
?scanBamFlag
?summarizeOverlaps

Parharn · 06-17-2014, 03:03 PM

Yes I got help in another thread, you are right. Thanks a lot though.
Here is the other thread just if any one wanted to know more.
https://www.biostars.org/p/103482/#103584

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
varScan output question	karl_s	Bioinformatics	0	01-09-2014 01:30 PM
the question from pairToBed output	g781	Bioinformatics	0	10-28-2013 12:51 AM
fastQC output question	arrchi	Bioinformatics	2	08-30-2012 02:06 PM
annovar output - question	dkrtndhkd	Bioinformatics	3	04-29-2012 10:23 AM
Question on BWA output	rboettcher	Bioinformatics	3	04-18-2011 04:21 AM

06-16-2014, 04:46 AM	#1
Parharn Member Location: Sweden Join Date: Jul 2013 Posts: 84	Zero output for sel.rn=rowSums(cnts) != 0 I get a problem when I want to try step 2 from "RNA-seq data pathway and gene-set analysis workflow" by Weijun Luo. I get zero rows with value !=0. Code: > p.cnts=assay(gnCnt) > cnts=p.cnts > dim(cnts) [1] 7017 2 > sel.rn=rowSums(cnts) !=0 > cnts=cnts[sel.rn,] > dim(cnts) [1] 0 2 When I check my BAM files with IGV it looks that I have mapped reads all over the genome. Any advise?

06-17-2014, 11:30 AM	#2
bigmw Senior Member Location: US Join Date: Aug 2013 Posts: 123	I haven’t seen your data and exactly how you did your reads mapping and counting. But please make sure on two things: First, you use the same version of reference genome at the TopHat step and in your TranscriptDb object. Second, you specify the right parameters in the reads counting step. They are different for single- vs paired-end data. The workflow example is for paired-end data. Single end data are treated differently: … flag <- scanBamFlag(isNotPrimaryRead=FALSE, isProperPair=NA) param <- ScanBamParam(flag=flag) gnCnt <- summarizeOverlaps(exByGn, bamfls, mode="Union", ignore.strand=TRUE, single.end=TRUE, param=param) … Please check help info for details: ?scanBamFlag ?summarizeOverlaps

06-17-2014, 03:03 PM	#3
Parharn Member Location: Sweden Join Date: Jul 2013 Posts: 84	Yes I got help in another thread, you are right. Thanks a lot though. Here is the other thread just if any one wanted to know more. https://www.biostars.org/p/103482/#103584