果蠅探索PRC復合物之文章復現(xiàn)（一）
老大曾在公眾號上發(fā)過一篇推文：在果蠅探索PRC復合物（逆向收費讀文獻2019-18），里面提到了
表觀領域的一篇文章：Global changes of H3K27me3 domains and Polycomb group protein distribution in the absence of recruiters Spps or Pho荸镊。這篇文章?lián)f是RNA-seq和ChIP-seq數(shù)據(jù)分析結(jié)合的典范料皇。關于PRC復合物,文章中的結(jié)論是-PcG recruiters, the PRC2 component E(z), and the PRC1 components Psc and Ph cobind thousands of active genes outside of H3K27me3 domains.關于如何獲取數(shù)據(jù)也就是拿到peaks文件和表達矩陣蝇裤，可以去嗶哩嗶哩網(wǎng)站搜索生信技能樹Chip—seq測序數(shù)據(jù)分析視頻。其實呢老大已經(jīng)錄制好了關于這篇文章中圖表復現(xiàn)的視頻葬荷，我這里跟著視頻喂窟，復現(xiàn)文章中的圖圖。
- 第一個復現(xiàn)的圖-Figure S12 Decreased expression of Pho and Spps in the corresponding mutants
上圖來自做完RNA-seq后的counts值椒舵，如果還不會linux的小伙伴又想直接練習R中的代碼蚂踊，可以問我要counts矩陣。
先要讀取表達矩陣笔宿，表達矩陣來自RNA-seq的feature counts的結(jié)果犁钟。也就是文件all.counts.id.txt。
讀進來的表達矩陣如下圖
- 第一列：基因名
- 第二列：染色體
- 第三列：基因起始坐標
- 第四列：基因終止坐標
- 第五列：略
為什么有些基因有這么多的其實坐標和終止坐標呢泼橘？是因為有些基因有多個外顯子涝动。
在復現(xiàn)上面的第一張圖片前，老大出題兒讓我對這個表達矩陣做進一步了解炬灭。
題目一
挑選出基因名前兩位為AB或者ab的基因名醋粟，并挑選出長度最長的基因
ab <- a[grep('^AB',a$Geneid,ignore.case = T),]
lapply(1:nrow(ab),function(i){
sum(as.numeric(strsplit(ab[i,4],';')[[1]])-as.numeric(strsplit(ab[i,3],';')[[1]]))
})
which.max(lapply(1:nrow(ab),function(i){
sum(as.numeric(strsplit(ab[i,4],';')[[1]])-as.numeric(strsplit(ab[i,3],';')[[1]]))
}))
ab[which.max(lapply(1:nrow(ab),function(i){
sum(as.numeric(strsplit(ab[i,4],';')[[1]])-as.numeric(strsplit(ab[i,3],';')[[1]]))
})),1]
或者可以
ab <- a[grep('^AB',a$Geneid,ignore.case = T),]
tmp2 <- apply(ab,1, function(i){
sum(as.numeric(strsplit(i[4],';')[[1]])-as.numeric(strsplit(i[3],';')[[1]]))
})
names(tmp) <- ab$Geneid#這一步是給向量tmp加上基因名，用names函數(shù)重归，給向量加基因名用names函數(shù)
tmp
names(which.max(tmp))
題目二：chr這一列有多少種元素
tmp3 <- as.data.frame(a$Chr)
unique(apply(tmp3,1, function(i){
unique(sort(str_split(i,';',simplify = T)[1,]))
}))
上面的代碼有點長米愿，需要分解一下
apply(tmp3,1, function(i){
unique(sort(str_split(i,';',simplify = T)[1,]))
})
得到的結(jié)果有17714行，實際上這個a矩陣就是有這么多鼻吮。那么我們在用uniq去重復就好了
現(xiàn)在對第一張圖復現(xiàn)育苟，圖是敲了Pho后看看敲除后的效果，rna-seq的表達量是否降低椎木。
rm(list = ls())
options(stringsAsFactors = F)
a=read.table('all.counts.id.txt',header = T)
dim(a)
cg1=a[a[,1]=='pho',7:16]
library(ggpubr)
library(stringr)
dat=data.frame(gene=as.numeric(cg),
sample=names(cg),
group=str_split(names(cg),'_',simplify = T)[,1]
)
ggbarplot(dat,x='sample',y='gene',color = 'group')
可以看到雖然上面的圖和文章的圖圖還有些差別违柏，關于作圖可以再另一個作圖的板塊細細的寫一下笨腥！上圖可以看出Pho敲除組的表達量是明顯低于其他組的，就OK了勇垛。
上面的代碼均來自生信技能樹哦
