提取
分别用rsem,kallisto,Salmon定量得到每个样本isoform的定量文件,提取了transcript_id一列和count一列,并且增加了序号列,使得每个transcript_id都有唯一的序号与之对应.以rsem定量结果中的Y43.isoforms.results为例
cd ~/rnaseq/rsem_out/ && cut -f 1,5 ~/rnaseq/rsem_out/Y43.isoforms.results >Y43_r && sed -i '1d' Y43_r && awk '$0=NR"\t"$0' Y43_r >Y43_r.csv && mv Y43_r.csv ~/rnaseq/compare/
cut -f 1,5 ~/rnaseq/rsem_out/Y43.isoforms.results >Y43_r：把Y43.isoforms.results的第1、5两列提取出来放到文件Y43_r中
写成bash文件,提取所有的样本.
提取结果
定量文件散点图
样本:Y43
red:rsem
blue:Salmon
green:kallisto
rsem&&Salmon
rsem&&kallisto
rsem&&kallisto
样本:Y45
rsem&&Salmon
rsem&&kallisto
Salmon&&kallisto
样本:O70
rsem&&Salmon
rsem&&kallisto
Salmon&&kallisto
样本:O77
rsem&&Salmon
rsem&&kallisto
Salmon&&kallisto
差异分析结果比对
rsem
Kallisto
Salmon
rsem
kallisto
Salmon
data <- read.table("Y48_s.csv") #读取文件
head(data)
查看
Y48_s =data$V3 #提取第三列
Y48_s=Y48_s+0.000001 #加上一个非常小的数防止后面做除法时出现分母为0的情况
data <- read.table("Y48_k.csv") #对由kallisto定量得到的结果采取一样的处理
Y48_k =data$V3
Y48_k=Y48_k+0.000001
data <- read.table("Y48_r.csv")
Y48_r =data$V3
Y48_r=Y48_r+0.000001
SK=Y48_s/ Y48_k
SR=Y48_s*2/ Y48_r
KR=Y48_k*3/ Y48_r
X=data$V1
par(pin=c(3,3))
plot(x=X, y=log10(SK), type="p", xlab="serial number",ylab="salmon/kallisto", main="SK",pch=20,cex=0.01)
plot(x=X, y=log10(SR), type="p", xlab="serial number",ylab="salmon/rsem", main="SR",pch=20,cex=0.01)
plot(x=X, y=log10(KR), type="p", xlab="serial number",ylab="kallisto/rsem", main="KR",pch=20,cex=0.01)
丑图
丑图
上面两种画图方式都不直观也不美观，数据量庞大，可能需要去除一些极端的值再看看数据的分布范围，做出合理的比较
以样本Y48为例
- 读入rsem和salmon的counts
- 用summary查看整体情况，初步了解
初步查看
-
去掉未表达的基因的影响
data<-data[data$V3!=0,]test<-test[test$V3!=0,]
删除0表达 -
去掉极端值
找出极端值
找极端 -
去掉极端值影响
去掉影响
1.931
1.68898488
1.638436
这个是两次处理前后的数据集均值比，data$mean/test$mean
去掉0值以后二者的统计结果明显更接近了，而最大值虽然相差很大，但是面对三万多行的数据，它的影响被弱化了。
下面通过画密度曲线观察定量结果的分布
初步尝试：将salmon的定量结果分成5个区间，目的是展示不同数量级的count的分布情况
-
困难：由于数据量过大，最大值高达几十万，而平均值仅仅200+或者100+，坐标轴的均匀划分，带来的问题是低值高分布区被挤压得与Y轴几近重合，高值低分布区则贴近X轴，如下：
看不到分布 -
为了解决这个问题，不再将一款软件的不同区间放在一起比较，而是将一个区间的不同软件放在一起比较
于是得到五个区间的不同软件定量结果比较图
0-10
分开展示
10-100
100-1000
1000-10000
10000—
这样比较貌似比散点图美观一些。
