哈爾濱醫科大學藥物轉錄組學實驗

為便於後輩使用，特意上傳！！

題目

1.利用T檢驗結合FC的方法識別差異表達蛋白編碼基因及lncRNA（ESC與CM之間的差異）注：考慮資料的預處理（例如 刪除全部樣本低表達基因，log轉化）。（結果檔案，基因ID ESC均值 CM均值 FC P值 FDR值）
2.繪製差異表達蛋白編碼基因及lncRNA火山圖（ggplot2）。
3.利用差異表達mRNA的表達譜進行聚類分析（pheatmap）
4.計算與差異表達lncRNA共表達的蛋白編碼基因（同時用斯皮爾曼相關及皮爾森相關），比較兩種方法得到的顯著共表達蛋白編碼基因列表的差異（cor.test）。
5.使用通過皮爾森計算得到的lncRNA（選擇差異表達程度最大的lncRNA）共表達蛋白編碼基因進行功能富集，從而預測lncRNA功能。
6.選擇差異表達程度最大的lncRNA，針對他的共表達蛋白編碼基因進行KEGG通路註釋。
擴充套件題
使用DEseq 計算蛋白編碼基因差異表達基因（利用read 數進行差異分析）

第一題

FLDR.IN <- "E:/1_生信相關/徐娟師姐/心血管發育相關_實驗課/ESC_CM"
FLDR.OUT <- "E:/1_生信相關/徐娟師姐/心血管發育相關_實驗課/result"


## 讀取mRNA的fpkm檔案
pro <- read.table(paste0(FLDR.IN, "/", "protein_coding.fpkm.landscape.txt"), sep = "\t", header = T, stringsAsFactors = F)
## 保留所有樣本中rpkm均大於1的基因
tmp <- apply(pro, 1, function(x){all(x > 1)})
protein <- pro[names(tmp[tmp == T]), ]  ## 7140 gene
## t.test
protein2 <- t(apply(protein, 1, function(x){log2(x+0.05)}))

pvalues <- NULL
for (i in 1:nrow(protein2)){
  pvalue <- t.test(protein2[i, grep("ESC", colnames(protein2))], protein2[i, grep("CM", colnames(protein2))])$p.value
  pvalues <- c(pvalues, pvalue)
}

fdr <- p.adjust(pvalues)

## fold change
mean_esc <- apply(protein[,grep("ESC",colnames(protein))], 1, mean)
mean_cm <- apply(protein[,grep("CM",colnames(protein))], 1, mean)
fc <- mean_cm / mean_esc
log2fc <- log2(fc)
## 合併所有結果
result <- data.frame(pvalues = pvalues, fdr = fdr, log2fc = log2fc, stringsAsFactors = F)
## 卡FDR小於0.01且fold change 大於2或者小於0.5
sig_result <- result[which(result$fdr < 0.01 & (result$log2fc > 1 | result$log2fc < -1)), ]  ## 3377 rows
## 加一列標籤
sig_result$significant <- "NULL"
sig_result[which(sig_result$log2fc > 1),4] <- "up"
sig_result[which(sig_result$log2fc < 1),4] <- "down"

## 輸出結果
write.table(sig_result, paste0(FLDR.OUT, "/", "result_ttest.txt"), sep = "\t", quote = F, col.names = T, row.names = T)

火山圖；熱圖

FLDR.IN <- "E:/1_生信相關/徐娟師姐/心血管發育相關_實驗課/ESC_CM"
FLDR.OUT <- "E:/1_生信相關/徐娟師姐/心血管發育相關_實驗課/result"
## 輸入結果檔案
result <- read.table(paste0(FLDR.OUT, "/", "result_ttest.txt"), sep = "\t", header = T, stringsAsFactors = F)
head(result)

## 火山圖 --------------------------------------------------
install.packages("ggplot2")
library(ggplot2)
library(ggplot2)
## 作圖
r03 <- ggplot(result, aes(log2fc, -log10(fdr)))
r03 + geom_point()
## 改變點的顏色
r03 + geom_point(color ="red")
r03 + geom_point(aes(color ="red"))
r03 + geom_point(aes(color = significant))

## 設定標題
r03xy <- r03 + geom_point(aes(color =significant))
r03xy + labs(title="Volcanoplot")
## 自定義顏色
r03xyp <- r03xy + labs(title="Volcanoplot")
r03xyp + scale_color_manual(values =c("#4876FF","#FF3E96"))

## 熱圖 -----------------------------------------------------
all_rpkm <- read.table(paste0(FLDR.IN, "/", "protein_coding.fpkm.landscape.txt"), sep = "\t", header = T, stringsAsFactors = F)
sig_rpkm <- all_rpkm[rownames(result), ]

## 定義標籤顏色
annotation_col <- data.frame(stage = factor(c(rep("ESC", 21), rep("CM", 16))))
rownames(annotation_col) <- colnames(sig_rpkm)
ann_colors <- list(stage = c(ESC = "#9BCD9B", CM = "#FFB5C5"))
## heatmap
ht <- as.matrix(sig_rpkm)
pheatmap(ht, show_colnames= T, show_rownames= F, scale= "row", fondsize= 6.5,
         annotation_col= annotation_col, 
         annotation_colors= ann_colors, 
         col = colorRampPalette(c("dodgerblue","white","red3"), alpha = T)(100),
         main = "heatmap of significant gene"
)

其他

data=read.csv("GPCR_DEG.txt",stringsAsFactors=F,sep="\t",skip=0,header =T)
data$threshold <- as.factor(ifelse(data$fdr < 0.05 & abs(log2(data$foldchange)) >=1,
                                   ifelse(log2(data$foldchange) > 1,'Up','Down'),'Not'))
library(ggplot2)
pdf("GPCR_DEG.pdf", height=6, width=8)
ggplot(data=data, 
       aes(x=log2(foldchange), y =-log10(pvalues), 
           colour=threshold,fill=threshold)) +
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_point(alpha=0.4, size=1.2) +
  xlim(c(-4, 4)) +
  theme_bw(base_size = 12, base_family = "Times") +
  geom_vline(xintercept=c(-1,1),lty=4,col="grey",lwd=0.6)+
  geom_hline(yintercept = -log10(0.05),lty=4,col="grey",lwd=0.6)+
  theme(legend.position="right",
        panel.grid=element_blank(),
        legend.title = element_blank(),
        legend.text= element_text(face="bold", color="black",family = "Times", size=8),
        plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(face="bold", color="black", size=12),
        axis.text.y = element_text(face="bold",  color="black", size=12),
        axis.title.x = element_text(face="bold", color="black", size=12),
        axis.title.y = element_text(face="bold",color="black", size=12))+
  labs(x="log2 (fold change)",y="-log10 (p-value)",title="GPCRs")
dev.off()

插入連結與圖片

圖片:

![Alt]
OVER!
其他的我也沒有了！！！！！！！！！！