2020.07.01-07.15学习小结

一、TCGA数据下载

1、TCGAbiolinks包

• 介绍：可用于检索，下载，并准备TCGA数据用于下游分析的R包
• 优点：具备一体化的下载整合，无需再使用复杂的方法对下载的单个数据重新进行整合。
• 使用

  GDCquery_clinic(project , type = "clinical") //临床数据下载
  qurey <- GDCquery(project, data.category, data.type, workflow.type,
  		legacy = FALSE, access, platform, file.type, barcode,
   		 experimental.strategy, sample.type) 	
  GDCdownload(query, method = "api", files.per.chunk = 100)
  data <- GDCprepare(query)	 
  
  	1. project
  		 TCGAbiolinks:::getGDCprojects()$project_id	//查看有哪些癌症类型
  	2. data.category
  		TCGAbiolinks:::getProjectSummary(project)	//查看某癌症有哪些数据类型
  	3. data.type
  		data.type = "Gene Expression Quantification"//下载rna-seq的counts数据
  		data.type = "miRNA Expression Quantification"//下载miRNA数据
  		data.type = "Copy Number Segment"			//下载Copy Number Variation数据
  		......
  	4. workflow.type
  		workflow.type="HTSeq - Counts"		//原始count数
  		workflow.type="HTSeq - FPKM-UQ"		//FPKM上四分位数标准化值
  		workflow.type="HTSeq - FPKM"		//FPKM值/表达量值
  		......
  	 5. legacy：不同的数据来源 Legacy 与 harmonized
  		  GDC Legacy Archive:以前在CGHUB和TCGA数据门户中存储的数据的原始数据，由TCGA数据协调中心（DCC）托管，
  	  在该门户中用GRCH37（HG19）和GRCH36（HG18）作为参考基因组	
  	   	  GDC harmonized database:可用数据与grch38（hg38）使用gdc生物信息学流程进行协调，该流程提供了生物标
  	  本和临床数据标准化的方法，简单讲就是对数据进行了一定标准化处理。harmonized数据库包括转录谱数据，甲基化数
  	  据，miRNA数据，但缺少芯片数据
  	6. 其他参数
  	  根据需求选择，可通过TCGA官网查看详细内容
  

二、免疫浸润

1、immunedeconv包

输入要求：行名为基因名(HGNC symbols )，列名为样本名的矩阵
使用：

res = deconvolute(gene_expression_matrix, method)

特殊情况

• CIBERSORT：

  //需要注明CIBERSORT.R与LM22.txt文件路径
   set_cibersort_binary("/path/to/CIBERSORT.R")
   set_cibersort_mat("/path/to/LM22.txt")
  

• TIMER:

  deconvolute(your_mixture_matrix, "timer",
  			indications=c("SKCM", "SKCM", "BLCA"))
  //indications是一个表示每个样本癌症类型的向量，TIMER支持的癌症类型查看可
  immunedeconv::timer_available_cancers
  

• xCell
  在cell_type_map中xCell 仅有39种细胞类型，故使用xCell网站运行：https://xcell.ucsf.edu/

三、生存分析

区别

• 统计描述
  估计生存率（生存函数）：Kaplan-Meier法

• 统计推断
  比较各组的生存率（生存曲线）：log-rank检验——单因素/分层分析
  生存期长短的影响因素分析： COX回归——单因素/多因素分析

K-M曲线可以很直观地表现出两组或多组的生存率或死亡率，因为生存数据往往都是非正态分布，因此常使用非参数的检验方法log-rank检验。在实际写作中，K-M曲线都会搭配log-rank检验的P值。

Cox回归本质上是一种回归模型，它没有直接使用生存时间，而是使用了危险度（hazard）作为因变量，该模型不用于估计生存率，而是用于因素分析，也就是找到某一个危险因素对结局事件发生的贡献度。

因此可使用K-M方法绘制生存曲线，使用log-rank检验评估生存差异。使用Cox比例风险模型进行多因素分析，评估预后因素作用。

1. 创建生存对象：Surv()
   详见 https://www.rdocumentation.org/packages/survival/versions/3.2-3/topics/Surv

   //useage
   Surv(time, time2, event,
       type=c('right', 'left', 'interval', 'counting', 'interval2', 'mstate'),
       origin=0)
   //example
   Surv(OSTime,OS)
   

2. 拟合生存曲线：survfit()
   详见https://www.rdocumentation.org/packages/survival/versions/3.2-3/topics/survfit

   //usegae
   survfit(formula, ...)
   //example
   survfit(Surv(OSTime, OS) ~ variable, data = data)
   

3. 绘制生存曲线：ggsurvplot()
   详见https://www.rdocumentation.org/packages/survminer/versions/0.4.7/topics/ggsurvplot

   //useage
   ggsurvplot(
     fit,					//survfit()对象/对象列表/包含生存曲线摘要的数据框
     data = NULL,			//用于拟合生存曲线的数据集，如未提供则从fit中提取
     fun = NULL,
     color = NULL,			//默认情况下，使用参数color =“ strata”按层为生存曲线着色，也可以使用参数调色板指定自定义调色板
     palette = NULL,		//调色板
     linetype = 1,			//线型
     conf.int = FALSE,		//绘制置信区间
     pval = FALSE,			//绘制P值
     pval.method = FALSE,	//Pval为T时，显示P值检验方法
     test.for.trend = FALSE,surv.median.line = "none",
     risk.table = FALSE,cumevents = FALSE,cumcensor = FALSE,
     tables.height = 0.25,group.by = NULL,facet.by = NULL,add.all = FALSE,
     combine = FALSE,ggtheme = theme_survminer(),tables.theme = ggtheme,
     ...
   )
   

4. 多因素分析：coxph()
   详见：https://www.rdocumentation.org/packages/survival/versions/3.2-3/topics/coxph

   coxph(formula,  //~左侧surv函数返回的生存对象，右侧为待分析因素
   	  data=, 	//用于解释formula/weights/subset中命名变量的数据框
   	  weights, subset, 
         na.action, init, control, 
         ties=c("efron","breslow","exact"), //默认为efron
         singular.ok=TRUE, robust, 
         model=FALSE, x=FALSE, y=TRUE, tt, method=ties,
         id, cluster, istate, statedata, ...)
   

四、 差异分析

1、使用FPKM数据进行差异基因分析

//将FPKM数据转换为TPM数据
expMatrix <- FPKM_matrix
fpkmToTpm <- function(fpkm)
{
  exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}
 apply(expMatrix,2,fpkmToTpm)

//限定分类顺序
Group <- factor(group,levels = ,ordered = F)
design2 <- model.matrix(~factor( Group ))

//表达矩阵数据校正
exprSet <-data
boxplot(exprSet,outline=FALSE, notch=T,col=Group, las=2)//查看数据是否需要标准化
library(limma) 
exprSet=normalizeBetweenArrays(exprSet)//数据标准化
boxplot(exprSet,outline=FALSE, notch=T,col=Group, las=2)
exprSet <- log2(exprSet+1)//需判断数据是否要log化

//差异分析：
dat <- exprSet
fit=lmFit(dat,design2)
fit=eBayes(fit)
options(digits = 4)
deg=topTable(fit,coef=2,adjust='BH',number = Inf)
difflimma <- subset(deg,abs(deg$logFC)>1&deg$P.Value<0.05)  //筛选条件为|logFC|>1，Pvalue<0.05

五、富集分析

library(clusterProfiler) 
library(org.Hs.eg.db)
//id转换
gene_id =bitr(gene_id,fromType = "SYMBOL",toType = "ENTREZID", OrgDb="org.Hs.eg.db")

1、GO分析

详见https://www.rdocumentation.org/packages/clusterProfiler/versions/3.0.4/topics/enrichGO

enrichGO(gene,    			    //entrez id
		 OrgDb, 				//组织数据库
		 keytype = "ENTREZID", 
		 ont = "MF", 			//MF", "BP", "CC"，"ALL"
		 pvalueCutoff = 0.05, 
		 pAdjustMethod = "BH", 
         universe, 
         qvalueCutoff = 0.2, 
         minGSSize = 10, 
         maxGSSize = 500, 
         readable = FALSE)

2、KEGG分析

详见https://www.rdocumentation.org/packages/clusterProfiler/versions/3.0.4/topics/enrichKEGG

enrichKEGG(gene, 
		  organism = "hsa", 
		  keyType = "kegg", 			//"kegg", 'ncbi-geneid', 'ncib-proteinid','uniprot'
		  pvalueCutoff = 0.05, 
		  pAdjustMethod = "BH", 
		  universe, 
		  minGSSize = 10, 
		  maxGSSize = 500, 
		  qvalueCutoff = 0.2, 
		  use_internal_data = FALSE)

六、ggplot

元素

1. plot：整张图，包括background和title
2. axis：包括stick，text，title
3. legend：包括background，text，title
4. facet：可分为外部stript（包括backgroud和text）和内部panel部分（包括backgroud、boder和网格线grid，其中粗的叫grid.major，细的叫grid.minor）。

函数类型

1. 用于运算（如fortify_，mean_等）
2. 初始化、展示（ggplot，plot，print等）
3. 按变量组图（facet_等）
4. 绘图（stat_，geom_，annotate），核心函数。
5. 微调图型：
   - scale_：与aes内的各种美学（shape、color、fill、alpha）调整有关的函数。
   - guides：调整所有的text。
   - coord_：调整坐标。
   - theme：调整不与数据有关的图的元素的函数。

1、ggpubr

介绍：基于ggplot2的用于绘制符合出版物要求的图形的可视化包
链接：https://www.jianshu.com/p/678213d605a5

分布图(Distribution)

#构建数据集
set.seed(1234)
df <- data.frame( sex=factor(rep(c("f", "M"), each=200)), 
weight=c(rnorm(200, 55), rnorm(200, 58)))
head(df)
##   sex   weight
## 1  f   53.79293
## 2  f   55.27743
## 3  f   56.08444
## 4  f   52.65430
## 5  f   55.42912
## 6  f   55.50606

密度分布图以及边际地毯线并添加平均值线

ggdensity(df, x="weight", add = "mean", rug = TRUE, color = "sex", fill = "sex",
palette = c("#00AFBB", "#E7B800"))

带有均值线和边际地毯线的直方图

gghistogram(df, x="weight", add = "mean", rug = TRUE, color = "sex", fill = "sex",
palette = c("#00AFBB", "#E7B800"))

箱线图与小提琴图

#加载数据集ToothGrowth
data("ToothGrowth")
df1 <- ToothGrowth
head(df1)
##    len  supp  dose
## 1  4.2   VC    0.5
## 2  11.5  VC    0.5
## 3  7.3   VC    0.5
## 4  5.8   VC    0.5
## 5  6.4   VC    0.5
## 6  10.0  VC    0.5
p <- ggboxplot(df1, x="dose", y="len", color = "dose", 
palette = c("#00AFBB", "#E7B800", "#FC4E07"), 
add = "jitter", shape="dose")#增加了jitter点，点shape由dose映射p

增加不同组间的p-value值，可以自定义需要标注的组间比较

my_comparisons <- list(c("0.5", "1"), c("1", "2"), c("0.5", "2"))
p+stat_compare_means(comparisons = my_comparisons)+#不同组间的比较
stat_compare_means(label.y = 50)

内有箱线图的小提琴图

ggviolin(df1, x="dose", y="len", fill = "dose", 
palette = c("#00AFBB", "#E7B800", "#FC4E07"), 
add = "boxplot", add.params = list(fill="white"))+ 
stat_compare_means(comparisons = my_comparisons, label = "p.signif")+#label这里表示选择显著性标记（星号） 
stat_compare_means(label.y = 50)

条形图

data("mtcars")
df2 <- mtcars
df2$cyl <- factor(df2$cyl)
df2$name <- rownames(df2)#添加一行name
head(df2[, c("name", "wt", "mpg", "cyl")])

按从小到大顺序绘制条形图

不分组排序
ggbarplot(df2, x="name", y="mpg", fill = "cyl", color = "white", 
palette = "jco",#杂志jco的配色 
sort.val = "desc",#下降排序 
sort.by.groups=FALSE,#不按组排序 
x.text.angle=60)
按组进行排序
ggbarplot(df2, x="name", y="mpg", fill = "cyl", color = "white", 
palette = "jco",#杂志jco的配色 
sort.val = "asc",#上升排序,区别于desc，具体看图演示 
sort.by.groups=TRUE,#按组排序 
x.text.angle=90)

偏差图
偏差图展示了与参考值之间的偏差

df2$mpg_z <- (df2$mpg-mean(df2$mpg))/sd(df2$mpg)
df2$mpg_grp <- factor(ifelse(df2$mpg_z<0, "low", "high"), levels = c("low", "high"))
head(df2[, c("name", "wt", "mpg", "mpg_grp", "cyl")])

绘制排序过的条形图

ggbarplot(df2, x="name", y="mpg_z", fill = "mpg_grp", color = "white", 
palette = "jco", sort.val = "asc", sort.by.groups = FALSE, x.text.angle=60, 
ylab = "MPG z-score", xlab = FALSE, legend.title="MPG Group")

坐标轴变换

ggbarplot(df2, x="name", y="mpg_z", fill = "mpg_grp", color = "white", 
palette = "jco", sort.val = "desc", sort.by.groups = FALSE, 
x.text.angle=90, ylab = "MPG z-score", xlab = FALSE, 
legend.title="MPG Group", rotate=TRUE, ggtheme = theme_minimal())

点图(Dot charts)
棒棒糖图(Lollipop chart)
棒棒图可以代替条形图展示数据

ggdotchart(df2, x="name", y="mpg", color = "cyl", 
palette = c("#00AFBB", "#E7B800", "#FC4E07"), sorting = "ascending", 
add = "segments", ggtheme = theme_pubr())

可以自设置各种参数

ggdotchart(df2, x="name", y="mpg", color = "cyl", 
palette = c("#00AFBB", "#E7B800", "#FC4E07"), sorting = "descending", 
add = "segments", rotate = TRUE, group = "cyl", dot.size = 6, 
label = round(df2$mpg), font.label = list(color="white", size=9, vjust=0.5), 
ggtheme = theme_pubr())

偏差图

ggdotchart(df2, x="name", y="mpg_z", color = "cyl", 
palette = c("#00AFBB", "#E7B800", "#FC4E07"), sorting = "descending", 
add = "segment", add.params = list(color="lightgray", size=2), 
group = "cyl", dot.size = 6, label = round(df2$mpg_z, 1), 
font.label = list(color="white", size=9, vjust=0.5), ggtheme = theme_pubr())+ 
geom_line(yintercept=0, linetype=2, color="lightgray")

Cleveland点图

ggdotchart(df2, x="name", y="mpg", color = "cyl", 
palette = c("#00AFBB", "#E7B800", "#FC4E07"), sorting = "descending", 
rotate = TRUE, dot.size = 2, y.text.col=TRUE, ggtheme = theme_pubr())+ 
theme_cleveland()

2、色板扩充

1. 默认颜色

   brewer.pal.info		#查询默认颜色
   

• seq类型：单渐变色，一种主色由浅到深
• qual类型：区分色，几种区分度很高的颜色组合
• div类型：双渐变色，一种颜色到另外一种颜色的渐变，有两种主色

2. 色板扩充

   //直接使用默认色板
   ggplot(...)+scale_fill_brewer(palette="Set1")+...
   
   //colorRamPalette使用默认色板
   library(RColorBrewer)
   colors <-colorRampPalette(brewer.pal(9,"Set1"))
   //colorRamPalette使用自定义色板
   colors <-colorRampPalette(c("BlueViolet","MediumOrchid","Orchid","MediumSlateBlue","CornflowerBlue",
   							"LightSteelBlue","LightSkyBlue","SkyBlue","PowDerBlue","PaleTurquoise",
   							"LightSeaGreen","Turquoise","LightGreen","DarkSeaGreen","Khaki",
   							"NavajoWhite","Gold","Orange"))
   ggplot(...+...+ scale_fill_manual(values = colors(22))+...)//22表示渐变梯度