吨吨吨的生信技能树r语言学习笔记（1）

前面课程的代码见
【生信技能树】R语言学习代码笔记 - 简书 (jianshu.com)

6.选取差异明显的基因的表达量矩阵绘制热图
感觉这一节课听的有点懵了

rm(list =ls())
options(stringAsFactors = F)
load(file = 'step1-output.Rdata')
dat=dat[,-3]
group_list=group_list[-3]
dat[1:4,1:4]
dim(dat)
dat=log2(dat)
dat[1:4,1:4]

d_h <- function(dat,group_list){
    cg=names(tail(sort(apply(dat,1,sd)),1000))##选取其中对dat所有行做标准差排序取最大的1000个基因的名称
    library(pheatmap)
    pheatmap(dat[cg,],show_colnames=F,show_rownames=F)#选取行名为cg的绘制热图
    n=t(scale(t(dat[cg,])))
    n[n>2]=2
    n[n<-2] = -2#标准化
    pheatmap(n,show_colnames=F,show_rownames=F)
    ac=data.frame(g=group_list)
    rownames(ac)=colnames(n)#分组
    pheatmap(n,show_colnames=F,show_rownames=F,
            annotation_col = ac)#绘图
}

d_h(dat,group_list)

dat <- dat[,-3]
group_list <- group_list[-3]

d_h(dat,group_list)

7. id转换
   为什么要进行id转换：关于基因ID的二三事 - 知乎 (zhihu.com)
   如何分离基因号
   法一： strsplit('分割内容'，‘分割符号’）应用循环找出
   法二：library（stringr)
   str_split(数据，'[.]' , simplify=T)变为matrix[,1]取列

library(org.Hs.eg.db)
导入a
g2s=toTable(org.Hs.egSYMBOL)
g2e=toTable(org.Hs.egENSEMBL)
b=merge(a,g2e, by='ensembl_id', all.x=T) #all.X=T,a所有的都给保留
d=merge(b,g2s, by='gene_id', all.x=T)
#现在开始处理NA
table(d$ensembl_id)[table(d$ensembl_id)>1]#查看重复序列
d=d[order(d$v1),] #调序
d=d[!duplicated(d$v1),]#去重
d=d[match(a$v1,d$v1),]#a中和d中的顺序相同
write.csv (d, 'file.csv')

8. 任意基因任意癌症表达量分组的生存分析
   介绍：
   TCGA数据库的初次了解 - 简书 (jianshu.com)
   (10条消息) GEO和TCGA_我是小飞熊的博客-CSDN博客_tcga和geo数据库的区别
   生存曲线＋单图复现，这个lncRNA数据库好用到飞起！还不快学！_腾讯新闻 (qq.com)
   手把手教你解读生存曲线 - 知乎 (zhihu.com)
   (10条消息) [转]生存分析之Kaplan-Meier曲线_为你千千万万遍z的博客-CSDN博客_kaplan-meier生存曲线
   www.oncolnc.org
   OncoLnc数据库收集了TCGA中21种肿瘤，共8647个病人的生存数据，以及对应的mRNA和miRNA的表达谱数据。同时收录了来自MiTranscriptome项目中lncRNA的表达谱数据，因此整个数据库包含了mRNA，miRNA，以及lncRNA共3种转录本的生存分析，方便用户挖掘各种肿瘤中和生存相关的基因。所以画个重点，别看这个数据库叫OncoLnc，事实上它的功能远不止lncRNA哦~

a <- read.csv("BRCA_7157_50_50.csv",sep = ",",header = T)
dat=a
library(ggplot2)
library(survival)
library(survminer) 
table(dat$Status)  #查看多少活得多少死的
dat$Status=ifelse(dat$Status=='Dead',1,0) #死了是1活着是0
sfit <- survfit(Surv(Days, Status)~Group, data=dat)#画生存曲线
sfit
summary(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)
ggsave('survival_TP53_in_BRCA_TCGA.png')

ggsurvplot(sfit,palette = c("#E7B800","#2E9FDF"),risk.table = T,
           pval = T,conf.int = T,xlab="time in months", 
           ggtheme =theme_light(),ncensor.plot=T )

#risk.table: TRUE or FALSE specifying whether to show or not the risk table. Default is FALSE.
#ncensor.plot:删失The height of the censor plot. Used when ncensor.plot = TRUE.
#pval:显示p值
#conf.int：logical value. If TRUE, plots confidence interval画置信区间.

画完了大概就是这个样子了，有一种不明觉厉的感觉

b119248752b60902d88c3b1a5c85f35.png

9.任意基因任意癌症表达量和临床性状关联
这节课使用了cbioportal
网站: http://www.cbioportal.org/index.do
介绍：TCGA数据挖掘终结者：cBioPortal - 知乎 (zhihu.com)

本次查看 Ovarian Serous Cystadenocarcinoma (TCGA, Nature 2011)中基因 ARHGAP18不同clinical stages表达是否有差异

a <- read.table("plot (1).txt",sep = '\t',header = T,fill = T)
colnames(a)=c('id','stage','gene','mut')
library(ggstatsplot)
ggbetweenstats(a,x = stage, y = gene)#琴型图
library(ggplot2)

可以得到结果

cafecc3f6ea15617b7dcd0777def3dc.png

10.表达矩阵的相关性

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")#下载airway数据包
BiocManager::install("airway") 
#edgeR包同理

what is airway?:
This package provides a Ranged Summarized Experiment object of read counts in genes for an RNA-Seq experiment on four human airway smooth muscle cell lines treated with dexamethasone. Details on the gene model and read counting procedure are provided in the package vignette.

library(airway)
data(airway)
exprSet <- assay(airway)  #基因的表达，列是处理情况行是基因
colnames(exprSet)
group_list=colData(airway)[,3]#提取处理情况分组
exprSet=exprSet[apply(exprSet,1,function(x)sum(x>1)>5),]
#提取exprSet的行中至少有五项大于0的部分
exprSet=log(edgeR::cpm(exprSet)+1) #计算差异取log+1
exprSet=exprSet[names(sort(apply(exprSet,1,mad),decreasing = T)[1:500]),]     
#将expeSet中行的中位数绝对偏差降序排序取前500个的基因名的表达信息
M=cor(log2(exprSet+1)) #计算相关性
tmp=data.frame(g=group_list)
rownames(tmp)=colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)#注释作图

结果

0663ff5836c1760c72fe22b411ffd01.png

11.芯片表达矩阵下游分析
几个应用的包 (limma的example)：
1.makeContrast:
Construct the contrast matrix corresponding to specified contrasts of a set of parameters。 (contrast matrix比较proges和stable)
2.lmFit:
Fit linear model for each gene given a series of arrays
3.contrasts.fit:
Given a linear model fit to microarray data, compute estimated coefficients and standard errors for a given set of contrasts.
4.eBayes
Empirical Bayes Statistics for Differential Expression
Given a linear model fit from lmFit, compute moderated t-statistics, moderated F-statistic, and log-odds of differential expression by empirical Bayes moderation of the standard errors towards a global value.差异表达结果

12.RNA-seq表达矩阵差异分析

rm(list = ls()) 
options(stringsAsFactors = F)
library(airway)
data("airway") 
exprSet=assay(airway)#获取表达矩阵
colnames(exprSet)
group_list=colData(airway)[,3]
exprSet=exprSet[apply(exprSet, 1,function(x)sum(x>1)>5),]
table(group_list)
exprSet[1:4,1:4]
boxplot(log(exprSet+1))