生信人的20個(gè)R語(yǔ)言習(xí)題 http://www.bio-info-trainee.com/3409.html
生信人的20個(gè)R語(yǔ)言習(xí)題 答案 http://www.bio-info-trainee.com/3415.html
參考答案：
https://cldiao.github.io/2018/09/27/R%E8%AF%AD%E8%A8%8020%E9%A2%98/
http://rvdsd.top/2018/08/30/R/20%E4%B8%AAR%E8%AF%AD%E8%A8%80%E4%B9%A0%E9%A2%98/
http://www.reibang.com/p/788010093c90

1. 安裝一些R包：

數(shù)據(jù)包： ALL, CLL, pasilla, airway
軟件包：limma，DESeq2，clusterProfiler
工具包：reshape2
繪圖包：ggplot2
不同領(lǐng)域的R包使用頻率不一樣髓考，在生物信息學(xué)領(lǐng)域，尤其需要掌握bioconductor系列包软能。

if(F){
source("http://bioconductor.org/biocLite.R")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")#修改鏡像沼沈，安裝會(huì)加速
BiocManager::install("clusterProfiler")
BiocManager::install("ComplexHeatmap")
BiocManager::install("maftools")
BiocManager::install("ggplot2")
BiocManager::install("jmzeng1314/biotrainee")
}
#或者如下：
source("https://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
BiocManager::install(c('ALL','CLL','pasilla','clusterProfiler')) 
BiocManager::install(c('airway','DESeq2','edgeR','limma'))
install.packages("reshape2", "ggplot2")

2.了解ExpressionSet對(duì)象

比如CLL包里面就有data(sCLLex) ，找到它包含的元素逸寓，提取其表達(dá)矩陣(使用exprs函數(shù))若治，查看其大小

1. 參考：http://www.bio-info-trainee.com/bioconductor_China/software/limma.html
2. 參考：https://github.com/bioconductor-china/basic/blob/master/ExpressionSet.md

suppressPackageStartupMessages(library(CLL))
data(sCLLex)
exprSet=exprs(sCLLex)
##sCLLex是依賴于CLL這個(gè)package的一個(gè)對(duì)象
samples=sampleNames(sCLLex)
pdata=pData(sCLLex)
group_list=as.character(pdata[,2])
dim(exprSet)
# [1] 12625    22
exprSet[1:5,1:5]
#            CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL
# 1000_at    5.743132  6.219412  5.523328  5.340477  5.229904
# 1001_at    2.285143  2.291229  2.287986  2.295313  2.662170
# 1002_f_at  3.309294  3.318466  3.354423  3.327130  3.365113
# 1003_s_at  1.085264  1.117288  1.084010  1.103217  1.074243
# 1004_at    7.544884  7.671801  7.474025  7.152482  6.902932

3.了解 str,head,help函數(shù)慨蓝，作用于 第二步提取到的表達(dá)矩陣

str(exprSet)
# str: Compactly display the internal structure of an R object, a diagnostic function and an alternative to summary (and to some extent, dput).
head(exprSet)

4. 安裝并了解 hgu95av2.db 包,看看 ls("package:hgu95av2.db") 后 顯示的那些變量

hgu95av2.db是一個(gè)注釋包，它為hgu95av2平臺(tái)的芯片提供注釋端幼，這個(gè)包中有很多注釋文件礼烈，如下所示：

BiocManager::install("hgu95av2.db")
library(hgu95av2.db)
ls("package:hgu95av2.db")
#[1] "hgu95av2"              "hgu95av2.db"           "hgu95av2_dbconn"       "hgu95av2_dbfile"       "hgu95av2_dbInfo"       "hgu95av2_dbschema"    
 [7] "hgu95av2ACCNUM"        "hgu95av2ALIAS2PROBE"   "hgu95av2CHR"           "hgu95av2CHRLENGTHS"    "hgu95av2CHRLOC"        "hgu95av2CHRLOCEND"    
[13] "hgu95av2ENSEMBL"       "hgu95av2ENSEMBL2PROBE" "hgu95av2ENTREZID"      "hgu95av2ENZYME"        "hgu95av2ENZYME2PROBE"  "hgu95av2GENENAME"     
[19] "hgu95av2GO"            "hgu95av2GO2ALLPROBES"  "hgu95av2GO2PROBE"      "hgu95av2MAP"           "hgu95av2MAPCOUNTS"     "hgu95av2OMIM"         
[25] "hgu95av2ORGANISM"      "hgu95av2ORGPKG"        "hgu95av2PATH"          "hgu95av2PATH2PROBE"    "hgu95av2PFAM"          "hgu95av2PMID"         
[31] "hgu95av2PMID2PROBE"    "hgu95av2PROSITE"       "hgu95av2REFSEQ"        "hgu95av2SYMBOL"        "hgu95av2UNIGENE"       "hgu95av2UNIPROT"      

5. 理解 head(toTable(hgu95av2SYMBOL)) 的用法，找到 TP53 基因?qū)?yīng)的探針I(yè)D

?hgu95av2SYMBOL
?toTable
summary(hgu95av2SYMBOL)
#SYMBOL map for chip hgu95av2 (object of class "ProbeAnnDbBimap")
|
| Lkeyname: probe_id (Ltablename: probes)
|    Lkeys: "1000_at", "1001_at", ... (total=12625/mapped=11460)
|
| Rkeyname: symbol (Rtablename: gene_info)
|    Rkeys: "A1BG", "A2M", ... (total=61050/mapped=8585)
|
| direction: L --> R

ids <- toTable(hgu95av2SYMBOL)
View(ids)
library(dplyr)
# 方法1：
filter(ids, symbol=="TP53") #用dplyr包的篩選功能婆跑，找到 TP53 基因?qū)?yīng)的探針I(yè)D
#     probe_id symbol
#1   1939_at   TP53
#2 1974_s_at   TP53
#3  31618_at   TP53

#方法2：
ids[grep("^TP53$", ids$symbol),]
#       probe_id symbol
# 966    1939_at   TP53
# 997  1974_s_at   TP53
# 1420  31618_at   TP53
#方法1此熬，2雖然結(jié)果相同，但是定義的對(duì)象是不同的

hug95av2SYMBOL是一個(gè)R對(duì)象，它提供的是芯片生產(chǎn)廠家與基因縮寫(xiě)之間的映射信息犀忱。這個(gè)映射的信息主要依據(jù)Entrez Gene數(shù)據(jù)庫(kù)∧蓟眩現(xiàn)在我們通過(guò)mappedkeys()這個(gè)函數(shù)，得到映射到基因上的探針信息阴汇。

6.理解探針與基因的對(duì)應(yīng)關(guān)系数冬，總共多少個(gè)基因，基因最多對(duì)應(yīng)多少個(gè)探針搀庶，是哪些基因拐纱，是不是因?yàn)檫@些基因很長(zhǎng)，所以在其上面設(shè)計(jì)多個(gè)探針呢哥倔？

length(unique(ids$symbol))
# [1] 8585
tail(sort(table(ids$symbol)))
#YME1L1  GAPDH INPP4A    MYB PTGER3  STAT1 
#     7      8      8      8      8      8 
table(sort(table(ids$symbol)))
#  1    2    3    4    5    6    7    8 
# 6555 1428  451  102   22   16    6    5 

image.png

不管是Agilent芯片秸架，還是Affymetrix芯片，上面設(shè)計(jì)的探針都非常短咆蒿。最長(zhǎng)的如Agilent芯片上的探針东抹，往往都是60bp，但是往往一個(gè)基因的長(zhǎng)度都好幾Kb沃测。因此一般多個(gè)探針對(duì)應(yīng)一個(gè)基因府阀，取最大表達(dá)值探針來(lái)作為基因的表達(dá)量。

7.第二步提取到的表達(dá)矩陣是12625個(gè)探針在22個(gè)樣本的表達(dá)量矩陣芽突，找到那些不在 hgu95av2.db 包收錄的對(duì)應(yīng)著SYMBOL的探針。

提示：有1165個(gè)探針是沒(méi)有對(duì)應(yīng)基因名字的董瞻。

%in% 邏輯判斷

用法 a %in% table
a值是否包含于table中寞蚌，為真輸出TURE，否者輸出FALSE

table(rownames(exprSet)) %in% ids$probe_id
# %in% is a more intuitive interface as a binary operator, which returns a logical vector indicating if there is a match or not for its left operand.
n_exprSet <- exprSet[!(rownames(exprSet) %in% ids$probe_id),]
dim(n_exprSet)
# [1] 1165   22
View(n_exprSet)
# These probes are not in the package.

8.過(guò)濾表達(dá)矩陣钠糊，刪除那1165個(gè)沒(méi)有對(duì)應(yīng)基因名字的探針挟秤。

方法1：%in% 邏輯判斷

exprSet <- exprSet[rownames(exprSet) %in% ids$probe_id, ]
dim(exprSet)
[1] 11460    22
View(exprSet)
# These probes are in the package.

方法2 mappedkeys() 映射關(guān)系

 length(hgu95av2SYMBOL)
[1] 12625
probe_map <- hgu95av2SYMBOL
length(probe_map)
[1] 12625
#全部的探針數(shù)目
# [1] 12625
probe_info <- mappedkeys(probe_map)
length(probe_info)
[1] 11460
#探針與基因產(chǎn)生映射的數(shù)目
gene_info <- as.list(probe_map[probe_info])
# 轉(zhuǎn)化為數(shù)據(jù)表
length(gene_info)
[1] 11460
gene_symbol <- toTable(probe_map[probe_info])
# 從hgu95av2SYMBOL文件中，取出有映射關(guān)系的探針抄伍，并生成數(shù)據(jù)框給gene_symbol
head(gene_symbol)
   probe_id  symbol
1   1000_at   MAPK3
2   1001_at    TIE1
3 1002_f_at CYP2C19
4 1003_s_at   CXCR5
5   1004_at   CXCR5
6   1005_at   DUSP1

mappedkeys用法示例艘刚，幫助理解。

library(hgu95av2.db)
  x <- hgu95av2GO
  x
  length(x)
  count.mappedkeys(x)
  x[1:3]
  links(x[1:3])

  ## Keep only the mapped keys
  keys(x) <- mappedkeys(x)
  length(x)
  count.mappedkeys(x)
  x # now it is a submap

  ## The above subsetting can also be achieved with
  x <- hgu95av2GO[mappedkeys(hgu95av2GO)]

  ## mappedkeys() and count.mappedkeys() also work with an environment
  ## or a list
  z <- list(k1=NA, k2=letters[1:4], k3="x")
  mappedkeys(z)
  count.mappedkeys(z)

  ## retrieve the set of primary keys for the ChipDb object named 'hgu95av2.db'
  keys <- keys(hgu95av2.db)
  head(keys)

9.整合表達(dá)矩陣截珍，多個(gè)探針對(duì)應(yīng)一個(gè)基因的情況下攀甚，只保留在所有樣本里面平均表達(dá)量最大的那個(gè)探針。

A. 提示岗喉，理解 tapply,by,aggregate,split 函數(shù) , 首先對(duì)每個(gè)基因找到最大表達(dá)量的探針秋度。
B. 然后根據(jù)得到探針去過(guò)濾原始表達(dá)矩陣

ids=ids[match(rownames(exprSet),ids$probe_id),]
head(ids)
exprSet[1:5,1:5]
tmp = by(exprSet,ids$symbol,function(x) rownames(x)[which.max(rowMeans(x))] )
probes = as.character(tmp)
exprSet=exprSet[rownames(exprSet) %in% probes ,]
dim(exprSet)
View(head(exprSet))

10.把過(guò)濾后的表達(dá)矩陣更改行名為基因的symbol，因?yàn)檫@個(gè)時(shí)候探針和基因是一對(duì)一關(guān)系了钱床。

rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]
exprSet[1:5,1:5]
#  CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL
# MAPK3    5.743132  6.219412  5.523328  5.340477  5.229904
# TIE1     2.285143  2.291229  2.287986  2.295313  2.662170
# CYP2C19  3.309294  3.318466  3.354423  3.327130  3.365113
# CXCR5    1.085264  1.117288  1.084010  1.103217  1.074243
# CXCR5    7.544884  7.671801  7.474025  7.152482  6.902932

library(reshape2)
exprSet_L=melt(exprSet)
colnames(exprSet_L)=c('probe','sample','value')
exprSet_L$group=rep(group_list,each=nrow(exprSet))
head(exprSet_L)
#  probe    sample    value    group
#1   MAPK3 CLL11.CEL 5.743132 progres.
#2    TIE1 CLL11.CEL 2.285143 progres.
#3 CYP2C19 CLL11.CEL 3.309294 progres.
#4   CXCR5 CLL11.CEL 1.085264 progres.
#5   CXCR5 CLL11.CEL 7.544884 progres.
#6   DUSP1 CLL11.CEL 5.083793 progres.
View(head(exprSet))

11. 對(duì)第10步得到的表達(dá)矩陣進(jìn)行探索荚斯，先畫(huà)第一個(gè)樣本的所有基因的表達(dá)量的boxplot,hist,density ， 然后畫(huà)所有樣本的 這些圖

1. 參考：http://bio-info-trainee.com/tmp/basic_visualization_for_expression_matrix.html
2. 理解ggplot2的繪圖語(yǔ)法，數(shù)據(jù)和圖形元素的映射關(guān)系

### ggplot2
library(ggplot2)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
print(p)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_violin()
print(p)
p=ggplot(exprSet_L,aes(value,fill=group))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4)
print(p)
p=ggplot(exprSet_L,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4)
print(p)
p=ggplot(exprSet_L,aes(value,col=group))+geom_density()
print(p)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
p=p+stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="red")
p=p+theme_set(theme_set(theme_bw(base_size=20)))
p=p+theme(text=element_text(face='bold'),axis.text.x=element_text(angle=30,hjust=1),axis.title=element_blank())
print(p)

image.png

image.png

image.png

image.png

image.png

image.png

12.理解統(tǒng)計(jì)學(xué)指標(biāo)mean,median,max,min,sd,var,mad并計(jì)算出每個(gè)基因在所有樣本的這些統(tǒng)計(jì)學(xué)指標(biāo)事期，最后按照mad值排序滥壕，取top 50 mad值的基因，得到列表兽泣。

注意：這個(gè)題目出的并不合規(guī)绎橘，請(qǐng)仔細(xì)看。

g_mean <- tail(sort(apply(exprSet,1,mean)),50)
g_median <- tail(sort(apply(exprSet,1,median)),50)
g_max <- tail(sort(apply(exprSet,1,max)),50)
g_min <- tail(sort(apply(exprSet,1,min)),50)
g_sd <- tail(sort(apply(exprSet,1,sd)),50)
g_var <- tail(sort(apply(exprSet,1,var)),50)
g_mad <- tail(sort(apply(exprSet,1,mad)),50)
g_mad
names(g_mad)

 [1] "DUSP5"   "IGFBP4"  "H1FX"    "ENPP2"   "FLNA"    "CLEC2B"  "TSPYL2"  "ZNF266"  "S100A9"  "NR4A2"   "TGFBI"   "ARF6"    "APBB2"   "VCAN"    "RBM38"  
[16] "CAPG"    "PLXNC1"  "RGS2"    "RNASE6"  "VAMP5"   "CYBB"    "GNLY"    "CCL3"    "OAS1"    "ENPP2"   "TRIB2"   "ZNF804A" "H1FX"    "IGH"     "JUND"   
[31] "SLC25A1" "PCDH9"   "VIPR1"   "COBLL1"  "GUSBP11" "S100A8"  "HBB"     "FOS"     "LHFPL2"  "FCN1"    "ZAP70"   "IGLC1"   "LGALS1"  "HBB"     "FOS"    
[46] "SLAMF1"  "TCF7"    "DMD"     "IGF2BP3" "FAM30A" 

13.根據(jù)第12步驟得到top 50 mad值的基因列表來(lái)取表達(dá)矩陣的子集撞叨，并且熱圖可視化子表達(dá)矩陣金踪。試試看其它5種熱圖的包的不同效果。

## heatmap
library(pheatmap)
choose_gene=names(tail(sort(apply(exprSet,1,mad)),50))
choose_matrix=exprSet[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))
pheatmap(choose_matrix)

image.png

14.取不同統(tǒng)計(jì)學(xué)指標(biāo)mean,median,max,mean,sd,var,mad的各top50基因列表牵敷，使用UpSetR包來(lái)看他們之間的overlap情況胡岔。

## UpSetR
# https://cran.r-project.org/web/packages/UpSetR/README.html
library(UpSetR)
g_all <- unique(c(names(g_mean),names(g_median),names(g_max),names(g_min),
names(g_sd),names(g_var),names(g_mad) ))
dat=data.frame(g_all=g_all,
g_mean=ifelse(g_all %in% names(g_mean) ,1,0),
g_median=ifelse(g_all %in% names(g_median) ,1,0),
g_max=ifelse(g_all %in% names(g_max) ,1,0),
g_min=ifelse(g_all %in% names(g_min) ,1,0),
g_sd=ifelse(g_all %in% names(g_sd) ,1,0),
g_var=ifelse(g_all %in% names(g_var) ,1,0),
g_mad=ifelse(g_all %in% names(g_mad) ,1,0)
)
upset(dat,nsets = 7)

image.png

15.在第二步的基礎(chǔ)上面提取CLL包里面的data(sCLLex) 數(shù)據(jù)對(duì)象的樣本的表型數(shù)據(jù)。

pdata=pData(sCLLex)
group_list=as.character(pdata[,2])
group_list
# [1] "progres." "stable"   "progres." "progres." "progres." "progres." "stable"   "stable"   "progres." "stable"   "progres." "stable"   "progres." "stable"  
# [15] "stable"   "progres." "progres." "progres." "progres." "progres." "progres." "stable"  
dim(exprSet)
# [1] 8585   22
exprSet[1:5,1:5]
#      CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL
MAPK3    5.743132  6.219412  5.523328  5.340477  5.229904
TIE1     2.285143  2.291229  2.287986  2.295313  2.662170
CYP2C19  3.309294  3.318466  3.354423  3.327130  3.365113
CXCR5    7.544884  7.671801  7.474025  7.152482  6.902932
DUSP1    5.083793  7.610593  7.631311  6.518594  5.059087

16.對(duì)所有樣本的表達(dá)矩陣進(jìn)行聚類并且繪圖枷餐，然后添加樣本的臨床表型數(shù)據(jù)信息(更改樣本名)

## hclust
colnames(exprSet)=paste(group_list,1:22,sep='')
# Define nodePar
nodePar <- list(lab.cex = 0.6, pch = c(NA, 19),
                cex = 0.7, col = "blue")
hc=hclust(dist(t(exprSet)))
par(mar=c(5,5,5,10))
plot(as.dendrogram(hc), nodePar = nodePar, horiz = TRUE)

image.png

17.對(duì)所有樣本的表達(dá)矩陣進(jìn)行PCA分析并且繪圖靶瘸，同樣要添加表型信息。

# install.packages("ggfortify")
library(ggfortify)
exprSet <- exprs(sCLLex)
df <- as.data.frame(t(exprSet))
df$group <- group_list
# autoplot uses ggplot2 to draw a particular plot for an object of a particular class in a single command.
autoplot(prcomp(df[,1:(ncol(df)-1)]), data=df, colour = 'group')

image.png

18.根據(jù)表達(dá)矩陣及樣本分組信息進(jìn)行批量T檢驗(yàn)毛肋，得到檢驗(yàn)結(jié)果表格

## t.test
dat = exprSet
group_list=as.factor(group_list)
group1 = which(group_list == levels(group_list)[1])
group2 = which(group_list == levels(group_list)[2])
dat1 = dat[, group1]
dat2 = dat[, group2]
dat = cbind(dat1, dat2)
pvals = apply(exprSet, 1, function(x){
  t.test(as.numeric(x)~group_list)$p.value
})
p.adj = p.adjust(pvals, method = "BH")
avg_1 = rowMeans(dat1)
avg_2 = rowMeans(dat2)
log2FC = avg_2-avg_1
DEG_t.test = cbind(avg_1, avg_2, log2FC, pvals, p.adj)
DEG_t.test=DEG_t.test[order(DEG_t.test[,4]),]
DEG_t.test=as.data.frame(DEG_t.test)
head(DEG_t.test)
#           avg_1    avg_2     log2FC        pvals     p.adj
36129_at 7.875615 8.791753  0.9161377 1.629755e-05 0.2057566
37676_at 6.622749 7.965007  1.3422581 4.058944e-05 0.2436177
33791_at 7.616197 5.786041 -1.8301554 6.965416e-05 0.2436177
39967_at 4.456446 2.152471 -2.3039752 8.993339e-05 0.2436177
34594_at 5.988866 7.058738  1.0698718 9.648226e-05 0.2436177
32198_at 4.157971 3.407405 -0.7505660 2.454557e-04 0.3516678

19.使用limma包對(duì)表達(dá)矩陣及樣本分組信息進(jìn)行差異分析怨咪，得到差異分析表格，重點(diǎn)看logFC和P值润匙，畫(huà)個(gè)火山圖(就是logFC和-log10(P值)的散點(diǎn)圖)诗眨。

# DEG by limma
suppressMessages(library(limma))
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast.matrix
##這個(gè)矩陣聲明，我們要把progres.組跟stable進(jìn)行差異分析比較
##step1
fit <- lmFit(exprSet,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix) ##這一步很重要孕讳，大家可以自行看看效果
fit2 <- eBayes(fit2) ## default no trend !!!
##eBayes() with trend=TRUE
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput)
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)

## volcano plot
DEG=nrDEG
logFC_cutoff <- with(DEG,mean(abs( logFC)) + 2*sd(abs( logFC)) )
DEG$change = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
                              ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
)
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
                    '\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',])
)

g = ggplot(data=DEG, aes(x=logFC, y=-log10(P.Value), color=change)) +
  geom_point(alpha=0.4, size=1.75) +
  theme_set(theme_set(theme_bw(base_size=20)))+
  xlab("log2 fold change") + ylab("-log10 p-value") +
  ggtitle( this_tile ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
  scale_colour_manual(values = c('blue','black','red')) ## corresponding to the levels(res$change)
print(g)

image.png

image.png

20.對(duì)T檢驗(yàn)結(jié)果的P值和limma包差異分析的P值畫(huà)散點(diǎn)圖匠楚，看看哪些基因相差很大。

### different P values
head(nrDEG)
head(DEG_t.test)
DEG_t.test=DEG_t.test[rownames(nrDEG),]
plot(DEG_t.test[,3],nrDEG[,1])
plot(DEG_t.test[,4],nrDEG[,4])
plot(-log10(DEG_t.test[,4]),-log10(nrDEG[,4]))

image.png

image.png

image.png

rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]
exprSet[1:5,1:5]
exprSet['GAPDH',]
exprSet['ACTB',]
exprSet['DLEU1',]
library(ggplot2)
library(ggpubr)
my_comparisons <- list(
  c("stable", "progres.")
)
dat=data.frame(group=group_list,
               sampleID= names(exprSet['DLEU1',]),
               values= as.numeric(exprSet['DLEU1',]))
ggboxplot(
  dat, x = "group", y = "values",
  color = "group",
  add = "jitter"
)+
  stat_compare_means(comparisons = my_comparisons, method = "t.test")

image.png

## heatmap
library(pheatmap)
choose_gene=head(rownames(nrDEG),25)
choose_matrix=exprSet[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))
pheatmap(choose_matrix)

image.png

生信技能樹(shù)公益視頻合輯：學(xué)習(xí)順序是linux厂财，r芋簿，軟件安裝，geo璃饱，小技巧与斤，ngs組學(xué)！
B站鏈接：https://m.bilibili.com/space/338686099
YouTube鏈接：https://m.youtube.com/channel/UC67sImqK7V8tSWHMG8azIVA/playlists
生信工程師入門(mén)最佳指南：https://mp.weixin.qq.com/s/vaX4ttaLIa19MefD86WfUA
學(xué)徒培養(yǎng)：https://mp.weixin.qq.com/s/3jw3_PgZXYd7FomxEMxFmw
生信技能樹(shù) - 簡(jiǎn)書(shū) http://www.reibang.com/u/d645f768d2d5