RNA-seq的数据往往需要剔除表达较低的基因,一般认为在半数以上样本中基因表达需要CPM>1,对于coding gene来说可以更为严格点,如CPM>2,而对于non-coding gene来说则可以适当降低点,如CPM>0.5。
###################04232019 发现更简洁用法
countdata <- read.table("yourdir/LIV.tsv",header = T,row.names = 1)
cpm <- t(t(countdata)/colSums(countdata) * 1000000)
row.names(cpm) <- row.names(countdata)
pass <- rowSums(cpm > 1)
pass <- pass[pass > (ncol(countdata) / 2)]
cpm <- as.data.frame(subset(cpm, rownames(cpm) %in% names(pass)))
###################
rt<-read.table("~/Desktop/OneDrive/GTEx/cpm.tsv",header = T,row.names = 1)#读取文档
MyFunction<-function(x){
ifelse(x>1,1,0)
}
#设置一个计数函数
keep_genes<-vector()#设置一个空的变量储存基因
row_names<-row.names(rt)
for (i in row_names){
cpm<-rt[i,]
count<-apply(cpm,2,MyFunction)
a<-sum(count)
if (a>=ncol(rt)/2)
keep_genes<-c(keep_genes,i)
}
rt_keep<-rt[keep_genes,]
write.csv(rt_keep,"cpm_1.csv")
其他标准也一样(如CPM>X),只要修改MyFunction中ifelse(x>X,1,0)。