哈尔滨医科大学药物转录组学实验

哈尔滨医科大学药物转录组学实验

为便于后辈使用，特意上传！！

题目

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!
其他的我也没有了！！！！！！！！！！