RNA-seq分析

rm(list = ls())
#清空环境变量
options(stringsAsFactors = F)
#字符串
a=read.table('all.id.txt',header = T)
#读入表达矩阵all.id.txt文件
tmp=a[1:14,1:7]
#取1-14行，1-7列
meta=a[,1:6]
#取所有行，1-6列
exprSet=a[,7:ncol(a)]
#从第七列开始取行，ncol(a)是说第七列之后有多少列，全部取出来
colnames(exprSet)
#看一下都有哪些样本
a2=exprSet[,'SRR1039516.hisat.bam']

library(airway)
data(airway)
exprSet=assay(airway)
colnames(exprSet)
a1=exprSet[,'SRR1039516']
group_list=colData(airway)[,3]

a2=data.frame(id=meta[,1],a2=a2)
a1=data.frame(id=names(a1),a1=as.numeric(a1))
#以下步骤用来去掉ENSEMBL.号后面的数字，这个问题我也遇到过啊
library(stringr)
a2$id <- str_split(a2$id,'\\.',simplify = T)[,1]

tmp=merge(a1,a2,by='id')
png('tmp.png')
plot(tmp[,2:3])
#只能看两个样本之间的相似性
dev.off()

library(corrplot)
png('cor.png')
corrplot(cor(log2(exprSet+1)))
#exprSet是表达矩阵，+1是为了防止log值为0；
dev.off()

library(pheatmap)
png('heatmap.png')
pheatmap(scale(cor(log2(exprSet+1))))
dev.off()