RNA-seq之QAPA量化后差异分析

RNA-seq之QAPA量化后差异分析

最近想看看AMD数据在转录层面有什么差异，所以先做了QAPA量化，从RNA-seq数据分析替代多聚腺苷酸化（APA）。 QAPA的安装和使用在GitHub上有详细的教程，附上链接 其实现在做APA量化分析的主流软件是 DaPars2，但是，我前面研究的时候，跑不通GitHub上的package，我真的急了什么都干得出来，在服务器上调试麻烦，我就拎到Windows上跑，还把人家一个.py文件拆开成两部分，在Windows上跑出中间文件temp_out，后面因为bedtools和它的Python“平替”——pybedtools都只能在linux上跑，所以用中间文件在服务器上跑了，估计除了我别人不会这样了！！但是废了这么大劲还是没戏，我觉得是我的数据下载的不对，DaPars2的Readme真的不如QAPA（ 附上链接），数据来源（仿照example在UCSC下载，但是没有属性名，只知道12列）。当然我用他的sample成功了，他用的hg19，现在大多数用hg38了。 之所以觉得是我的数据下载不对，是因为我在UCSC上又下载了hg19，跑数据的时候有error，和我在hg38同样的error，淦！大无语事件！搞不懂！

希望DaPars2配成功的宝贝可以私信教教我~

差异分析——求pvalue

1、数据集

给大家一个test数据，提取码：qapa
qapa quant 后的数据大概长这样（下文中这个文件：D:/RData/qapa_df/GSE99287/pau_results_99287.txt）

这么看真别扭，看下面这个吧

具体每一步就不写了，注释写在code里面，大家自己读一读吧！**看完觉得so easy哦！**我平时用Python多，但是做这个的时候觉得R更快，我几乎不用，都是现查现学，所以整体效率还欠佳，但是后续用的多了知道的多一些可能会有提升！这个就当是练手了。

2、思路流程

希望大家养成好的代码习惯，在解决问题之前先写好框架（本科老师以前教过伪代码，但是我全忘记了，不记得格式什么的，所以很粗糙的写个大概每一步要做的），这个可以帮助我缕清思路。任何时候拿到task先不着急上手，想逻辑，逻辑整清楚，知道每一步做什么，再上手，节省时间而且高效！

2.1 读数据

2.2 去除多余属性，得到拼接矩阵

拼接矩阵格式：
gene_name：normal_PAU：amd_PAU：normal_TPM：amd_TPM

2.3 卡TPM值

（1）在amd和normal中每个都至少存在3个TPM>1的值
（2）卡完TPM后将拼接矩阵分割，不要TPM列

2.4 秩和检验求pvalue

（1）秩和检验求p值
（2）求矫正后的p，即p.adjust
（3）卡pvalue<0.05，存为resCOX

2.5 t 检验求pvalue

（1）删除一行中全部相同的rows
（2）t检验求p值
（3）求矫正后的p，即p.adjust
（4）卡pvalue<0.05，存为resT

2.6 存储数据

重要结果直接存本地，下次直接看，不用重新浪费时间跑了！（懒人最大的福音！我就是）

3、Coding

setwd('D:/RData/qapa_df/GSE99287')
#读入qapa的结果
qapa<-read.table('D:/RData/qapa_df/GSE99287/pau_results.txt',header = T,row.names =1)#读文件
normal<-read.table('D:/RData/qapa_df/GSE99287/normal.txt')#读normal数据，正常数据作为对照组
amd<-read.table('D:/RData/qapa_df/GSE99287/amd.txt')#读amd数据
sample_num=nrow(normal)+nrow(amd)#所有样本数
#创建normal和amd的column-name
for(i in 1:nrow(normal)){
  normal[i,2]=paste(normal[i,1],"_quant.PAU",sep = "")
  normal[i,3]=paste(normal[i,1],"_quant.TPM",sep = "")
}
for(i in 1:nrow(amd)){
  amd[i,2]=paste(amd[i,1],"_quant.PAU",sep = "")
  amd[i,3]=paste(amd[i,1],"_quant.TPM",sep = "")
}
#将原数据提取出基因名，PAU和TPM拼接为新的数据
#得到的数据大概是gene_NAME：nature_PAU:amd_PAU:nature_TPM:amd_TPM
qdata<-qapa[3]#基因名
#normal_PAU
cols_normal=which(colnames(qapa) %in% normal[,2])
for(i in 1:length(cols_normal)){
  qdata<-cbind(qdata,qapa[cols_normal[i]])
}
#amd_PAU
cols_amd=which(colnames(qapa) %in% amd[,2])
for(i in 1:length(cols_amd)){
  qdata<-cbind(qdata,qapa[cols_amd[i]])
}
#normal_TPM
cols_normal=which(colnames(qapa) %in% normal[,3])
for(i in 1:length(cols_normal)){
  qdata<-cbind(qdata,qapa[cols_normal[i]])
}
#amd_TPM
cols_amd=which(colnames(qapa) %in% amd[,3])
for(i in 1:length(cols_amd)){
  qdata<-cbind(qdata,qapa[cols_amd[i]])
}
qdata<-na.omit(qdata)#删除有NA的所有行
#卡TPM值,TPM>1.normal 和 amd中都至少有三个样本的TPM>1
pv_data=qdata
delrow_tpm<-c()
j<-1
tpm_normal_start=sample_num+2
tpm_normal_end=tpm_normal_start+nrow(normal)-1
tpm_amd_start=tpm_normal_end+1
tpm_amd_end=ncol(pv_data)
for(i in 1:nrow(pv_data)){
  normal_tpm_count=sum(colSums(qdata[i,tpm_normal_start:tpm_normal_end]>1))
  amd_tpm_count=sum(colSums(qdata[i,tpm_amd_start:tpm_amd_end]>1))
  if(normal_tpm_count<3||amd_tpm_count<3){
    delrow_tpm[j]<-i
    j<-j+1
  }
}
pv_data<-pv_data[-delrow_tpm,]
sltedData<-pv_data[,1:(sample_num+1)]
#做检验，求p值
gene_num<-c(1:nrow(sltedData))
normal_start=2
normal_end=2+nrow(normal)-1
amd_start=normal_end+1
amd_end=ncol(sltedData)

p.w<-rep(NA,nrow(sltedData))
for(i in gene_num){
  p.w[i]<- wilcox.test(as.numeric(sltedData[i,normal_start:normal_end]),as.numeric(sltedData[i,amd_start:amd_end]),paired = FALSE)$p.value;
}
#秩和检验结果，卡p值
p.w<- as.numeric(p.w)
fdr.w<-p.adjust(p.w,method="fdr",length(p.w)) 
resCox<-cbind(sltedData[1],p.w,fdr.w)
resCox<-subset(resCox,p.w<0.05)

#t检验数据处理，全0或全100都删除
tdata=sltedData
tdata$sum=rowSums(tdata[,2:ncol(tdata)])
delcol<-c()
j<-1

for(i in gene_num){
  if(tdata[i,60]==0||tdata[i,60]==5800){
    delcol[j]<-i
    j<-j+1
  }
}
tdata<-tdata[-delcol,]
p.t<-rep(NA,nrow(tdata))
for(i in 1:nrow(tdata)){
  p.t[i]<- t.test(tdata[i,normal_start:normal_end],tdata[i,amd_start:amd_end],paired = FALSE)$p.value;
}
#t检验数据处理，拼接后，p-adjust
fdr.t<-p.adjust(p.t,method="fdr",length(p.t)) 
resT<-cbind(tdata[1],p.t,fdr.t)
resT <- subset(resT, p.t < 0.05 )

write.csv(resCox,file = 'different_gene.csv')

4、坑！巨坑！——字符串拼接（paste）

normal[i,2]=paste(normal[i,1],"_quant.PAU",sep = "")

paste(a,b)函数默认的拼接的时候会在a和b之间加空格！！一定要加sep=“”，才是"无缝衔接"

5、小tips

5.1 矩阵拼接

cbind()列拼接：左右
rbind()行拼接：上下

5.2 判断一行有几列满足特定条件

sum(colSums(qdata[5,1:20]>1))

第5行从第1列到第20列数据中大于1的列数。好用的！

结

886，去吃饭了，回来跑着editing数据再写一篇富集分析！