前面課程的代碼見
【生信技能樹】R語言學(xué)習(xí)代碼筆記 - 簡書 (jianshu.com)
6.選取差異明顯的基因的表達(dá)量矩陣?yán)L制熱圖
感覺這一節(jié)課聽的有點(diǎn)懵了
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))##選取其中對(duì)dat所有行做標(biāo)準(zhǔn)差排序取最大的1000個(gè)基因的名稱
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#標(biāo)準(zhǔn)化
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)
- id轉(zhuǎn)換
為什么要進(jìn)行id轉(zhuǎn)換：關(guān)于基因ID的二三事 - 知乎 (zhihu.com)
如何分離基因號(hào)
法一： strsplit('分割內(nèi)容'，‘分割符號(hào)’）應(yīng)用循環(huán)找出
法二：library（stringr)
str_split(數(shù)據(jù)，'[.]' , simplify=T)變?yōu)閙atrix[,1]取列
library(org.Hs.eg.db)
導(dǎo)入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)
#現(xiàn)在開始處理NA
table(d$ensembl_id)[table(d$ensembl_id)>1]#查看重復(fù)序列
d=d[order(d$v1),] #調(diào)序
d=d[!duplicated(d$v1),]#去重
d=d[match(a$v1,d$v1),]#a中和d中的順序相同
write.csv (d, 'file.csv')
- 任意基因任意癌癥表達(dá)量分組的生存分析
介紹：
TCGA數(shù)據(jù)庫的初次了解 - 簡書 (jianshu.com)
(10條消息) GEO和TCGA_我是小飛熊的博客-CSDN博客_tcga和geo數(shù)據(jù)庫的區(qū)別
生存曲線＋單圖復(fù)現(xiàn)，這個(gè)lncRNA數(shù)據(jù)庫好用到飛起！還不快學(xué)！_騰訊新聞 (qq.com)
手把手教你解讀生存曲線 - 知乎 (zhihu.com)
(10條消息) [轉(zhuǎn)]生存分析之Kaplan-Meier曲線_為你千千萬萬遍z的博客-CSDN博客_kaplan-meier生存曲線
www.oncolnc.org
OncoLnc數(shù)據(jù)庫收集了TCGA中21種腫瘤，共8647個(gè)病人的生存數(shù)據(jù)，以及對(duì)應(yīng)的mRNA和miRNA的表達(dá)譜數(shù)據(jù)。同時(shí)收錄了來自MiTranscriptome項(xiàng)目中l(wèi)ncRNA的表達(dá)譜數(shù)據(jù)，因此整個(gè)數(shù)據(jù)庫包含了mRNA，miRNA，以及l(fā)ncRNA共3種轉(zhuǎn)錄本的生存分析，方便用戶挖掘各種腫瘤中和生存相關(guān)的基因。所以畫個(gè)重點(diǎn)，別看這個(gè)數(shù)據(jù)庫叫OncoLnc，事實(shí)上它的功能遠(yuǎn)不止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畫置信區(qū)間.
畫完了大概就是這個(gè)樣子了，有一種不明覺厲的感覺
9.任意基因任意癌癥表達(dá)量和臨床性狀關(guān)聯(lián)
這節(jié)課使用了cbioportal
網(wǎng)站: http://www.cbioportal.org/index.do
介紹：TCGA數(shù)據(jù)挖掘終結(jié)者：cBioPortal - 知乎 (zhihu.com)
本次查看 Ovarian Serous Cystadenocarcinoma (TCGA, Nature 2011)中基因 ARHGAP18不同clinical stages表達(dá)是否有差異
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)
可以得到結(jié)果
10.表達(dá)矩陣的相關(guān)性
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")#下載airway數(shù)據(jù)包
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) #基因的表達(dá)，列是處理情況行是基因
colnames(exprSet)
group_list=colData(airway)[,3]#提取處理情況分組
exprSet=exprSet[apply(exprSet,1,function(x)sum(x>1)>5),]
#提取exprSet的行中至少有五項(xiàng)大于0的部分
exprSet=log(edgeR::cpm(exprSet)+1) #計(jì)算差異取log+1
exprSet=exprSet[names(sort(apply(exprSet,1,mad),decreasing = T)[1:500]),]
#將expeSet中行的中位數(shù)絕對(duì)偏差降序排序取前500個(gè)的基因名的表達(dá)信息
M=cor(log2(exprSet+1)) #計(jì)算相關(guān)性
tmp=data.frame(g=group_list)
rownames(tmp)=colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)#注釋作圖
結(jié)果
11.芯片表達(dá)矩陣下游分析
幾個(gè)應(yīng)用的包 (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.差異表達(dá)結(jié)果
12.RNA-seq表達(dá)矩陣差異分析
rm(list = ls())
options(stringsAsFactors = F)
library(airway)
data("airway")
exprSet=assay(airway)#獲取表達(dá)矩陣
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))
