Cellranger使用教程

建库,人和小鼠的数据库可以直接下载,对于无法直接下载的需要自行下载全基因组序列和gtf文件,根据 cellranger mkref构建参考数据库

wget ftp://ftp.ensembl.org/pub/release-97/fasta/danio_rerio/dna/Danio_rerio.GRCz11.dna.primary_assembly.fa.gz
gunzip Danio_rerio.GRCz11.dna.primary_assembly.fa.gz
wget ftp://ftp.ensembl.org/pub/release-97/gtf/danio_rerio/Danio_rerio.GRCz11.97.gtf.gz
gunzip Danio_rerio.GRCz11.97.gtf.gz
cellranger mkgtf Danio_rerio.GRCz11.97.gtf Danio_rerio.GRCz11.97.filtered.gtf                  
                 --attribute=gene_biotype:protein_coding \
                 --attribute=gene_biotype:lincRNA \
                 --attribute=gene_biotype:antisense \
                 --attribute=gene_biotype:IG_LV_gene \
                 --attribute=gene_biotype:IG_V_gene \
                 --attribute=gene_biotype:IG_V_pseudogene \
                 --attribute=gene_biotype:IG_D_gene \
                 --attribute=gene_biotype:IG_J_gene \
                 --attribute=gene_biotype:IG_J_pseudogene \
                 --attribute=gene_biotype:IG_C_gene \
                 --attribute=gene_biotype:IG_C_pseudogene \
                 --attribute=gene_biotype:TR_V_gene \
                 --attribute=gene_biotype:TR_V_pseudogene \
                 --attribute=gene_biotype:TR_D_gene \
                 --attribute=gene_biotype:TR_J_gene \
                 --attribute=gene_biotype:TR_J_pseudogene \
                 --attribute=gene_biotype:TR_C_gene
cellranger mkref --nthreads=80 --genome=ref_zebr_GRCz11 --fasta=Danio_rerio.GRCz11.dna.primary_assembly.fa --genes=Danio_rerio.GRCz11.97.filtered.gtf --ref-version=3.1.0    #注意gene_id和transcript_id的顺序

计算表达量,下机的原始数据可以通过bcl2fastq拆分,也可通过cellranger自带的mkfastq命令拆分。bcl2fastq拆分后得到两个fastq文件,index信息包含在第一个fastq文件第一行尾部;mkfastq拆分后得到三个fastq文件,index信息包含在一个单独的fastq文件里。两种拆分结果都可作为count命令的输入,但文件的命名一定要严格安装软件的说明,否则会出错。

cellranger count --id=PBLzebr1 --transcriptome=/home/wuchangsong/sc_cell/ref_zebr_GRCz11 --fastqs=/home/wuchangsong/sc_cell/fastq/ --sample=PBLzebr1

得到的表达矩阵在filtered_feature_bc_matrix文件夹中,在analysis文件夹中有PCA和tSNE聚类结果文本形式。cellranger count 计算的结果只能作为错略观测的结果,如果需要进一步分析聚类细胞,还需要进行下游分析,这里使用官方推荐 R 包(Seurat),后边的分析参考Seurat的使用。

参考:https://support.10xgenomics.com/single-cell-gene-expression/software/release-notes/build#GRCh38mm10_3.1.0

Seurat使用流程

seurat软件安装

Depends R (>= 3.4.0), methods

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("Seurat")

CentOS系统安装时要注意gcc的版本

setwd("D:/Experiment_data/zxj/outs")
library(Seurat)
pbl.data <- Read10X(data.dir = "D:/Experiment_data/zxj/outs/filtered_feature_bc_matrix")
dim(pbl.data) #查看行和列
#创建 Seurat 对象与数据过滤。保留在>=3 个细胞中表达的基因;保留能检测到>=200 个基因的细胞。
pbl <- CreateSeuratObject(counts = pbl.data, project = "pbl1907", min.cells = 3, min.features = 200) 
#mt-开头的为线粒体基因,这里将其进行标记并统计其分布频率
pbl[["percent.mt"]] <- PercentageFeatureSet(pbl, pattern = "^mt-")
# 对 pbmc 对象做小提琴图,分别为基因数,细胞数和线粒体占比
VlnPlot(object = pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
#根据图片中基因数和线粒体数,分别设置过滤参数,这里基因数 200-2500,线粒体百分比为小于 5%
pbl <- subset(pbl, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")+ NoLegend()
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")+ NoLegend()
CombinePlots(plots = list(plot1, plot2))
#标准化
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
#鉴定高变基因
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
plot3 <- VariableFeaturePlot(pbmc)+ NoLegend()
plot4 <- LabelPoints(plot = plot3, points = top10, repel = TRUE,xnudge=0,ynudge=0)
CombinePlots(plots = list(plot1, plot2))
#数据归一化,对所有基因进行标准化,默认只是标准化高变基因( 2000 个),速度更快,不影响 PCA 和分群,但影响热图的绘制。
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc,features = all.genes,vars.to.regress = "percent.mt")
#线性降维(PCA),默认用高变基因集,但也可通过 features 参数自己指定
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#定义可视化细胞和功能的几种有用的方式PCA,包括VizDimReduction,DimPlot,和DimHeatmap
VizDimLoadings(object = pbmc, dims = 1:2, reduction = "pca")
DimPlot(pbmc, reduction = "pca")+ NoLegend()
DimHeatmap(pbmc, dims = 1:2, cells = 500, balanced = TRUE)
#鉴定数据集的可用维度,主成分分析结束后需要确定哪些主成分所代表的基因可以进入下游分析,这里可以使用JackStraw做重抽样分析。可以用JackStrawPlot可视化看看哪些主成分可以进行下游分析。
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot(pbmc, dims = 1:15)#虚线以上的为可用维度,你也可以调整 dims 参数,画出所有 pca 查看
#肘部图(碎石图),基于每个主成分对方差解释率的排名。建议尝试选择多个主成分个数做下游分析,对整体影响不大;在选择此参数时,建议选择偏高的数字( “宁滥勿缺”,为了获取更多的稀有分群);有些亚群很罕见,如果没有先验知识,很难将这种大小的数据集与背景噪声区分开来。
ElbowPlot(object = pbmc)
#非线性降维( UMAP/tSNE)
#基于 PCA 空间中的欧氏距离计算 nearest neighbor graph,优化任意两个细胞间的距离权重(输入上一步得到的 PC 维数)
pbmc <- FindNeighbors(pbmc, dims = 1:10)
#resolution 参数决定下游聚类分析得到的分群数,对于 3K 左右的细胞,设为 0.4-1.2 能得到较好的结果(官方说明);如果数据量增大,该参数也应该适当增大;增加的值会导致更多的群集。
pbmc <- FindClusters(pbmc, resolution = 0.5)
#使用 Idents()函数可查看不同细胞的分群;
head(Idents(pbmc), 8)
#Seurat 提供了几种非线性降维的方法进行数据可视化(在低维空间把相似的细胞聚在一起),比如 UMAP 和 t-SNE,运行 UMAP 需要先安装'umap-learn'包,两种方法都可以使用,但不要混用,这样,后面的结算结果会将先前的聚类覆盖掉,只能保留一个
# If you haven't installed UMAP, you can do so via reticulate::py_install(packages = 'umap-learn')
pbmc <- RunTSNE(pbmc, dims = 1:10)
#pbmc <- RunUMAP(object = pbmc, dims = 1:10)
#DimPlot(object = pbmc, reduction = "umap")
#用 DimPlot()函数绘制散点图, reduction = "tsne",指定绘制类型;如果不指定,默认先从搜索 umap, 然后 tsne, 再然后 pca;也可以直接使用这 3 个函数 PCAPlot()、 TSNEPlot()、UMAPPlot(); cols, pt.size 分别调整分组颜色和点的大小;
tsneplot<-TSNEPlot(pbmc,label = TRUE, pt.size = 1.5)+ NoLegend()
#保存数据
saveRDS(pbmc, file = "pbmc_tutorial.rds")
#寻找差异表达的特征(聚类生物标志物)
#寻找某个聚类和其他所有聚类显著表达的基因
cluster1.markers <- FindMarkers(object = pbmc, ident.1 = 1, min.pct = 0.25)
head(x = cluster1.markers, n = 5)
#find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(object = pbmc, ident.1 = 5, ident.2 = c(0,3), min.pct = 0.25)
# find markers for every cluster compared to all remaining cells, report only the positive ones
#pbmc.markers <- FindAllMarkers(object = pbmc, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25)
pbmc.markers <- FindAllMarkers(object = pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25, return.thresh = 0.01)
library(dplyr)
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(object = pbmc, features = top10$gene) + NoLegend()
#画图
VlnPlot(object = pbmc, features = c("MS4A1", "CD79A"))
FeaturePlot(object = pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
# you can plot raw UMI counts as well
VlnPlot(object = pbmc, features.plot = c("NKG7", "PF4"), use.raw = TRUE, y.log = TRUE)
FeaturePlot(object = pbmc, features.plot = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"), cols.use = c("grey", "blue"), reduction.use = "tsne")
#将单元类型标识分配给集群
new.cluster.ids <- c("Memory CD4 T", "CD14+ Mono", "Naive CD4 T", "B", "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Mk")
names(x = new.cluster.ids) <- levels(x = pbmc)
pbmc <- RenameIdents(object = pbmc, new.cluster.ids)
tsneplot2<-TSNEPlot(pbmc,label = TRUE, pt.size = 1.5)+ NoLegend()





未完待续。。。

awk命令将fastq文件转换成fasta

awk '{if(NR%4==1||NR%4==2){print $0}}'  test.fq | sed 's/^@/>/g' > myfile.fasta

 有时候,比如说我们需要做一个clustalw,我们需要将多行的fasta文件(multiple line fasta)格式文件,转换成单行的fasta文件(single line fasta)序列文件,怎么办呢。 巧用awk + printf一行搞定。
awk '/^>/ { print n $0;}  !/^>/ {printf "%s", $0, n="\n"}  END {print ""}'  test.fa

 

WGCNA(my project)

setwd("D:/Ma_transcription/WGCNA/")
database <- read.table(file = "COS_deseq_counts_normalized.txt", sep = "\t", header = T, row.names = 1, stringsAsFactors = F)
library(reshape2)
library('WGCNA')
enableWGCNAThreads()#打开多线程
WGCNA_matrix = t(database[order(apply(database,1,mad), decreasing = T)[1:10000],])
subname=sapply(colnames(database),function(x) strsplit(x,"_")[[1]][1])
datTraits = data.frame(gsm=names(database),
subtype=subname)
rownames(datTraits)=datTraits[,1]
head(datTraits)
# Choose a set of soft-thresholding powers
powers = c(c(1:10), seq(from = 12, to=20, by=2))
datExpr <- WGCNA_matrix
# Call the network topology analysis function
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
# Plot the results:
par(mfrow = c(1,2));
cex1 = 0.9;
# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
     main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=cex1,col="red");
# this line corresponds to using an R^2 cut-off of h
abline(h=0.90,col="green")
# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
     main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
sft$powerEstimate
net = blockwiseModules(datExpr, power = sft$powerEstimate,
                       maxBlockSize = 5000,TOMType = "unsigned", 
                       minModuleSize = 30,reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       saveTOMFileBase = "AS-green-FPKM-TOM",
                       verbose = 3)
table(net$colors)

# Convert labels to colors for plotting
mergedColors = labels2colors(net$colors)
table(mergedColors)
# Plot the dendrogram and the module colors underneath
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
## assign all of the gene to their corresponding module 
## hclust for the genes.
#明确样本数和基因数
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
#首先针对样本做个系统聚类树
datExpr_tree<-hclust(dist(datExpr), method = "average")
par(mar = c(0,5,2,0))
plot(datExpr_tree, main = "Sample clustering", sub="", xlab="", cex.lab = 2, 
     cex.axis = 1, cex.main = 1,cex.lab=1)
## 如果这个时候样本是有性状,或者临床表型的,可以加进去看看是否聚类合理
#针对前面构造的样品矩阵添加对应颜色
sample_colors <- numbers2colors(as.numeric(factor(datTraits$subtype)), 
                                colors = c("grey","blue","red","green"),signed = FALSE)
## 这个给样品添加对应颜色的代码需要自行修改以适应自己的数据分析项目。
#  sample_colors <- numbers2colors( datTraits ,signed = FALSE)
## 如果样品有多种分类情况,而且 datTraits 里面都是分类信息,那么可以直接用上面代码,
##当然,这样给的颜色不明显,意义不大。
#构造10个样品的系统聚类树及性状热图
par(mar = c(1,4,3,1),cex=0.8)
plotDendroAndColors(datExpr_tree, sample_colors,
                    groupLabels = colnames(sample),
                    cex.dendroLabels = 0.8,
                    marAll = c(1, 4, 3, 1),
                    cex.rowText = 0.01,
design=model.matrix(~0+ datTraits$subtype)
colnames(design)=levels(datTraits$subtype)
moduleColors <- labels2colors(net$colors)
# Recalculate MEs with color labels
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0); ##不同颜色的模块的ME值矩阵(样本vs模块)
moduleTraitCor = cor(MEs, design , use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

sizeGrWindow(10,6)
# Will display correlations and their p-values
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                   signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = names(design),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = FALSE,
               colors = greenWhiteRed(50),
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 0.5,
               zlim = c(-1,1),
               main = paste("Module-trait relationships"))                 main = "Sample dendrogram and trait heatmap")
head(design)
# names (colors) of the modules
modNames = substring(names(MEs), 3)
geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"));
## 算出每个模块跟基因的皮尔森相关系数矩阵
## MEs是每个模块在每个样本里面的值
## datExpr是每个基因在每个样本的表达量
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));
names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");
## 只有连续型性状才能只有计算
## 这里把是否属于 CB 表型这个变量用0,1进行数值化。
CB = as.data.frame(design[,2]);
names(CB) = "CB"
geneTraitSignificance = as.data.frame(cor(datExpr, CB, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
names(geneTraitSignificance) = paste("GS.", names(CB), sep="");
names(GSPvalue) = paste("p.GS.", names(CB), sep="")
module = "turquoise"
column = match(module, modNames);
moduleGenes = moduleColors==module;
sizeGrWindow(7, 7);
par(mfrow = c(1,1));
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
                   abs(geneTraitSignificance[moduleGenes, 1]),
                   xlab = paste("Module Membership in", module, "module"),
                   ylab = "Gene significance for CB",
                   main = paste("Module membership vs. gene significance\n"),
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)

#首先针对所有基因画热图
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
geneTree = net$dendrograms[[1]]; 
#生成全基因不相似TOM矩阵
dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6); 
plotTOM = dissTOM^7; 
diag(plotTOM) = NA; 
#TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")

#然后随机选取部分基因作图
nSelect = 400
# For reproducibility, we set the random seed
set.seed(10);
select = sample(nGenes, size = nSelect);
selectTOM = dissTOM[select, select];
# There’s no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select];
# Open a graphical window
sizeGrWindow(9,9)
# Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing
# the color palette; setting the diagonal to NA also improves the clarity of the plot
plotDiss = selectTOM^7;
diag(plotDiss) = NA;
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")

#最后画模块和性状的关系
# Recalculate module eigengenes
MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
## 只有连续型性状才能只有计算
## 这里把是否属于 Luminal 表型这个变量用0,1进行数值化。
CB = as.data.frame(design[,2]);
names(CB) = "CB"
# Add the weight to existing module eigengenes
MET = orderMEs(cbind(MEs, CB))
# Plot the relationships among the eigengenes and the trait
sizeGrWindow(5,7.5);
par(cex = 0.9)
plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8,                         xLabelsAngle= 90)
# Plot the dendrogram
sizeGrWindow(6,6);
par(cex = 1.0)
## 模块的聚类图
plotEigengeneNetworks(MET, "Eigengene dendrogram", marDendro = c(0,4,2,0),
                      plotHeatmaps = FALSE)
# Plot the heatmap matrix (note: this plot will overwrite the dendrogram plot)
par(cex = 1.0)
## 性状与模块热图
plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2),
                      plotDendrograms = FALSE, xLabelsAngle = 90)
# Recalculate topological overlap
TOM = TOMsimilarityFromExpr(datExpr, power = 6); 
# Select module
module = "blue";
# Select module probes
probes = colnames(datExpr) ## 我们例子里面的probe就是基因名
inModule = (moduleColors==module);
modProbes = probes[inModule]; 
## 也是提取指定模块的基因名
# Select the corresponding Topological Overlap
modTOM = TOM[inModule, inModule];
dimnames(modTOM) = list(modProbes, modProbes)
nTop = 200;
IMConn = softConnectivity(datExpr[, modProbes]);
top = (rank(-IMConn) <= nTop)
vis = exportNetworkToCytoscape(modTOM[top, top], edgeFile = paste("CytoscapeInput-edges-top200", paste(module, collapse="-"), ".txt", sep=""), nodeFile = paste("CytoscapeInput-nodes-top200", paste(module, collapse="-"), ".txt", sep=""), weighted = TRUE, threshold = 0);

参考:lncRNA实战项目-第六步-WGCNA相关性分析

 

DESeq2使用流程

library("DESeq2")
database <- read.table(file = "COS_genes.count_table.xls", sep = "\t", header = T, row.names = 1)
countData <- database[,1:15]
condition <- factor(c(rep("COS0",3),rep("COS1",3),rep("COS3",3),rep("COS6",3),rep("COS12",3)), levels = c("COS0", "COS1","COS3","COS6","COS12"))
countData <- round(as.matrix(countData))
coldata <- data.frame(row.names = colnames(countData), condition)
dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), design= ~ condition )
dds <- dds[ rowSums(counts(dds)) > 1, ]
dds <- DESeq(dds)
res0vs1 <- results(dds, contrast = c("condition","COS1","COS0"))
table(res0vs1$pvalue <0.001)
res0vs1 <- res0vs1[order(res0vs1$pvalue),]
write.table(as.data.frame(res0vs1), file="res0vs1.xls", sep="\t", quote = F)

GSEA命令行

前边已经说到MSigDB数据库的entrez ID和symble ID都是人源的,其他物种要想做GSEA的话,就必须要ID转换,做成属于自己研究物种的数据库,这一步是关键。想要根据简单的Nr注释或其他注释来转换是行不通的,比如我研究的物种为硬骨鱼类,Nr注释结果分值最高的前20个一般也都是鱼类,就算注释到了人源的基因,也可能存在版本不对的问题。因此我们要知道MSigDB数据库基因ID的版本,然后下载对应的全基因组蛋白序列,然后比对卡阈值,做成自己的数据库。

GSEA的运行可以界面操作和命令行操作,界面操作教程很多就不多做赘述,它的不足之处除了操作繁琐之外,里面几个较大的库需要的内存较高,像我的4G渣渣电脑,好几个都运行失败,不得不转到服务器上运行;而且界面也不能批量运行。

java -cp gsea-3.0.jar -Xmx10000m xtools.gsea.Gsea -res ma/normalized_counts_all.gct -cls ma/samples.cls#COS1_versus_COS0 -gmx ma/ma.h.all.v6.1.symbols.gmt -collapse false -mode Max_probe -norm meandiv -nperm 1000 -permute gene_set -rnd_type no_balance -scoring_scheme weighted -rpt_label COS1vsCOS0_h -metric Ratio_of_Classes -sort real -order descending -create_gcts false -create_svgs false -include_only_symbols true -make_sets true -median false -num 100 -plot_top_x 20 -rnd_seed timestamp -save_rnd_lists false -set_max 2000 -set_min 5 -zip_report false -out multiple -gui false

-Xmx10000m:设置最大内存10000兆,之前界面操作失败就因为这个参数
-res:表达量文件
-cls:样本信息
-gmx:数据库文件
-nperm:Number of permutations
-rpt_label:输出文件的前缀
-metric:排序的方法,如果有重复,可以考虑使用T-test;无重复,可以考虑使用Ratio of calsses(差异倍数)或Diff of classes(差异绝对值)
-plot_top_x:默认20,代表富集分析排名最高的20个通路
-set_max:通路基因的最大数量,默认500,但由于某些通路基因数大于500,建议提高阈值
-set_min:通路基因的最小数量,默认15
-out:输出的文件夹
其他参数均为默认

版权声明:本文为博主原创文章,未经博主允许不得转载。

GSEA分析流程

什么是GSEA?

基因集富集分析 (Gene Set Enrichment Analysis, GSEA)用来评估一个预先定义的基因集的基因在与表型相关度排序的基因列表中的分布趋势,从而判断其对表型的贡献。基因集合富集分析检测基因集合而不是单个基因的表达变化,因此可以包含这些细微的表达变化,预期得到更为理想的结果。

GSEA原理

给定一个排序的基因表L和一个预先定义的基因集S (比如编码某个代谢通路的产物的基因, 基因组上物理位置相近的基因,或同一GO注释下的基因),GSEA的目的是判断S里面的成员s在L里面是随机分布还是主要聚集在L的顶部或底部。这些基因排序的依据是其在不同表型状态下的表达差异,若研究的基因集S的成员显著聚集在L的顶部或底部,则说明此基因集成员对表型的差异有贡献,也是我们关注的基因集。

GSEA分析中几个关键的概念

计算富集得分 (ES, enrichment score),ES反应基因集成员s在排序列表L的两端富集的程度。计算方式是,从基因集L的第一个基因开始,计算一个累计统计值。当遇到一个落在s里面的基因,则增加统计值。遇到一个不在s里面的基因,则降低统计值。每一步统计值增加或减少的幅度与基因的表达变化程度(更严格的是与基因和表型的关联度)是相关的。富集得分ES最后定义为最大的峰值。正值ES表示基因集在列表的顶部富集,负值ES表示基因集在列表的底部富集。

评估富集得分(ES)的显著性,通过基于表型而不改变基因之间关系的排列检验 (permutation test)计算观察到的富集得分(ES)出现的可能性。若样品量少,也可基于基因集做排列检验 (permutation test),计算p-value。

多重假设检验矫正,首先对每个基因子集s计算得到的ES根据基因集的大小进行标准化得到Normalized Enrichment Score (NES)。随后针对NES计算假阳性率。(计算NES也有另外一种方法,是计算出的ES除以排列检验得到的所有ES的平均值)

Leading-edge subset,对富集得分贡献最大的基因成员。

MSigDB数据库

MSigDB(Molecular Signatures Database)数据库中定义了已知的基因集合:http://software.broadinstitute.org/gsea/msigdb

包括H和C1-C7八个系列(Collection),每个系列内容为:H: hallmark gene sets(效应)特征基因集合,共50组(比如细胞凋亡特征基因集);C1: positional gene sets 位置基因集合,根据染色体位置,共326个;C2: curated gene sets:(专家)共识基因集合,基于通路、文献等(这部分包括我们熟悉的KEGG信号通路);C3: motif gene sets:模式基因集合,主要包括microRNA和转录因子靶基因两部分;C4: computational gene sets:计算基因集合,通过挖掘癌症相关芯片数据定义的基因集合;C5: GO gene sets:Gene Ontology 基因本体论,包括BP(生物学过程biological process,细胞原件cellular component和分子功能molecular function三部分);C6: oncogenic signatures:癌症特征基因集合,大部分来源于NCBI GEO 未发表芯片数据;C7: immunologic signatures: 免疫相关基因集合。

软件下载:http://software.broadinstitute.org/gsea/downloads.jsp

数据格式:http://www.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats

一:准备数据(GSEA 排序需要三个文件,如下)

1.表达量文件(*.gct)

2.说明文件(*.cls)

3. 基因集合(.gmt)

or

MSigDB提供的数据库是人源的,所以你的研究物种不是人的话,就需要自己制作gmt文件,这也是GSEA分析的难点和关键所在

DESeq分析流程

source("https://bioconductor.org/biocLite.R")
biocLite("DESeq")
library("DESeq")
database <- read.table(file = "macrophage_genes.count_table.matrix", sep = "\t", header = T, row.names = 1)
countData <- database[,1:3]
condition <- factor(c("A","B","C"))#无生物学重复
#type <- factor(c(rep("A",3), rep("B",3))) 有生物学重复
database <- round(as.matrix(countData))#取整数型
dds <- newCountDataSet(database,condition)
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds, method="blind", sharingMode="fit-only" )#无生物学重复
dds <- estimateDispersions(dds) # 有生物学重复
resAvsB <- nbinomTest(dds,"A","B")
resAvsC <- nbinomTest(dds,"A","C")
table(resAvsB$pval <0.05)
table(resAvsC$pval <0.05)
resAvsB <- resAvsB[order(resAvsB$pval),]#排序
resAvsC <- resAvsC[order(resAvsC$pval),] 
plotMA(resAvsB, ylim = c(-5,5), col = ifelse(resAvsB$pval>=0.05, "gray32", "red3"),linecol = "#ff000080" )#画MA图
norcounts <- counts(dds, normalized=T)#提取标准化后的counts
write.table(as.data.frame(norcounts),
file="normalized_counts_ABC.txt",
sep="\t",
quote = F)
resAvsBup <- subset(resAvsB, pval < 0.05 & log2FoldChange >0 )
resAvsBdown <- subset(resAvsC, pval < 0.05 & log2FoldChange <0 )

WGCNA分析流程

WGCNA是一个R包,对一个完全不会R的人来说,确实费了一番功夫,不过也将我对R的学习提前提上日程。

分析步骤:

1.WGCNA安装

source("http://bioconductor.org/biocLite.R")
biocLite(c("AnnotationDbi", "impute", "GO.db", "preprocessCore"))
install.packages("WGCNA")

2.输入数据准备

原始数据https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE48213,共56个样本;对数据的合并可参考http://www.biotrainee.com:8080/thread-603-1-1.html,输入文件格式如下:

 

setwd('WGCNA/')
fpkm=read.table('test.txt',sep = '\t',stringsAsFactors = F)
head(fpkm)
dim(fpkm) ##共56个细胞系,36953个基因的表达矩阵
names(fpkm)##细胞系名称

3.数据预处理

library(reshape2)
library('WGCNA')
RNAseq_voom <- fpkm
WGCNA_matrix = t(RNAseq_voom[order(apply(RNAseq_voom,1,mad), decreasing = T)[1:3000],]) ##top 3000 mad genes,t参数将矩阵变为符合WGCNA要求的形式:行名为gene,列名为样品
datExpr <- WGCNA_matrix
save(datExpr, file = "AS-green-FPKM-01-dataInput.RData")

4.确定最佳beta值

 

powers = c(c(1:10), seq(from = 12, to=20, by=2))
# Call the network topology analysis function
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
# Plot the results:
##sizeGrWindow(9, 5)
par(mfrow = c(1,2));
cex1 = 0.9;
# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
     main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=cex1,col="red");
# this line corresponds to using an R^2 cut-off of h
abline(h=0.90,col="red")
# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
     main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")

最佳的beta值就是sft$powerEstimate,此例为6

5.一步法构建共表达矩阵

net = blockwiseModules(datExpr, power = 6, maxBlockSize = 6000,
                       TOMType = "unsigned", minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       saveTOMFileBase = "AS-green-FPKM-TOM",
                       verbose = 3)
table(net$colors)
#绘图结果展示
mergedColors = labels2colors(net$colors)
# Plot the dendrogram and the module colors underneath
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
#结果保存
moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
table(moduleColors)
MEs = net$MEs;
geneTree = net$dendrograms[[1]];
save(MEs, moduleLabels, moduleColors, geneTree,
     file = "AS-green-FPKM-02-networkConstruction-auto.RData")

#对样本做个系统聚类树
#明确基因数和样本数
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
datExpr_tree<-hclust(dist(datExpr), method = "average")
par(mar = c(0,5,2,0))
plot(datExpr_tree, main = "Sample clustering", sub="", xlab="", cex.lab = 2, 
     cex.axis = 1, cex.main = 1,cex.lab=1)

6.导出网络到Cytoscape

# Recalculate topological overlap if needed
TOM = TOMsimilarityFromExpr(datExpr, power = 6);
# Read in the annotation file
# annot = read.csv(file = "GeneAnnotation.csv");
# Select modules需要修改,选择需要导出的模块颜色
modules = c("turquoise", "blue");
# Select module probes选择模块探测
probes = colnames(datExpr)
inModule = is.finite(match(moduleColors, modules));
modProbes = probes[inModule];
#modGenes = annot$gene_symbol[match(modProbes, annot$substanceBXH)];
# Select the corresponding Topological Overlap
modTOM = TOM[inModule, inModule];
dimnames(modTOM) = list(modProbes, modProbes)
# Export the network into edge and node list files Cytoscape can read
cyt = exportNetworkToCytoscape(modTOM,
                               edgeFile = paste("AS-green-FPKM-One-step-CytoscapeInput-edges-", paste(modules, collapse="-"), ".txt", sep=""),
                               nodeFile = paste("AS-green-FPKM-One-step-CytoscapeInput-nodes-", paste(modules, collapse="-"), ".txt", sep=""),
                               weighted = TRUE,
                               threshold = 0.02,
                               nodeNames = modProbes,
                               #altNodeNames = modGenes,
                               nodeAttr = moduleColors[inModule]);

7.分析网络可视化

用heatmap可视化权重网络,heatmap每一行或列对应一个基因,颜色越深表示有较高的邻近

options(stringsAsFactors = FALSE);
lnames = load(file = "AS-green-FPKM-01-dataInput.RData");
lnames
lnames = load(file = "AS-green-FPKM-02-networkConstruction-auto.RData");
lnames
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
#可视化全部基因网络
# Calculate topological overlap anew: this could be done more efficiently by saving the TOM
# calculated during module detection, but let us do it again here.
dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6);
# Transform dissTOM with a power to make moderately strong connections more visible in the heatmap
plotTOM = dissTOM^7;
# Set diagonal to NA for a nicer plot
diag(plotTOM) = NA;
# Call the plot function
#sizeGrWindow(9,9)
TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
#随便选取1000个基因来可视化
nSelect = 1000
# For reproducibility, we set the random seed
set.seed(10);
select = sample(nGenes, size = nSelect);
selectTOM = dissTOM[select, select];
# There's no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select];
# Open a graphical window
#sizeGrWindow(9,9)
# Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing
# the color palette; setting the diagonal to NA also improves the clarity of the plot
plotDiss = selectTOM^7;
diag(plotDiss) = NA;
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")

8.多步法构建网络

#2.选择合适的阀值,同上
#3. 网络构建:(1) Co-expression similarity and adjacency
softPower = 6;
adjacency = adjacency(datExpr, power = softPower);
#(2) 邻近矩阵到拓扑矩阵的转换,Turn adjacency into topological overlap
TOM = TOMsimilarity(adjacency);
dissTOM = 1-TOM
# (3) 聚类拓扑矩阵
#Call the hierarchical clustering function
geneTree = hclust(as.dist(dissTOM), method = "average");
# Plot the resulting clustering tree (dendrogram)
#sizeGrWindow(12,9)
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
     labels = FALSE, hang = 0.04);
#(4) 聚类分支的休整dynamicTreeCut
# We like large modules, so we set the minimum module size relatively high:
minModuleSize = 30;
# Module identification using dynamic tree cut:
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
                            deepSplit = 2, pamRespectsDendro = FALSE,
                            minClusterSize = minModuleSize);
table(dynamicMods)
#4. 绘画结果展示
# Convert numeric lables into colors
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)
# Plot the dendrogram and colors underneath
#sizeGrWindow(8,6)
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05,
                    main = "Gene dendrogram and module colors")
#5. 聚类结果相似模块的融合,Merging of modules whose expression profiles are very similar
#在聚类树中每一leaf是一个短线,代表一个基因,
#不同分之间靠的越近表示有高的共表达基因,将共表达极其相似的modules进行融合
# Calculate eigengenes
MEList = moduleEigengenes(datExpr, colors = dynamicColors)
MEs = MEList$eigengenes
# Calculate dissimilarity of module eigengenes
MEDiss = 1-cor(MEs);
# Cluster module eigengenes
METree = hclust(as.dist(MEDiss), method = "average");
# Plot the result
#sizeGrWindow(7, 6)
plot(METree, main = "Clustering of module eigengenes",
     xlab = "", sub = "")
#选择有75%相关性的进行融合
MEDissThres = 0.25
# Plot the cut line into the dendrogram
abline(h=MEDissThres, col = "red")
# Call an automatic merging function
merge = mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)
# The merged module colors
mergedColors = merge$colors;
# Eigengenes of the new merged modules:
mergedMEs = merge$newMEs;
#绘制融合前(Dynamic Tree Cut)和融合后(Merged dynamic)的聚类图
#sizeGrWindow(12, 9)
#pdf(file = "Plots/geneDendro-3.pdf", wi = 9, he = 6)
plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors),
                    c("Dynamic Tree Cut", "Merged dynamic"),
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
#dev.off()
# 只是绘制融合后聚类图
plotDendroAndColors(geneTree,mergedColors,"Merged dynamic",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
#5.结果保存
# Rename to moduleColors
moduleColors = mergedColors
# Construct numerical labels corresponding to the colors
colorOrder = c("grey", standardColors(50));
moduleLabels = match(moduleColors, colorOrder)-1;
MEs = mergedMEs;
# Save module colors and labels for use in subsequent parts
save(MEs, moduleLabels, moduleColors, geneTree, file = "AS-green-FPKM-02-networkConstruction-stepByStep.RData")
#6. 导出网络到Cytoscape
# Recalculate topological overlap if needed
TOM = TOMsimilarityFromExpr(datExpr, power = 6);
# Read in the annotation file
# annot = read.csv(file = "GeneAnnotation.csv");
# Select modules需要修改
modules = c("brown", "red");
# Select module probes
probes = colnames(datExpr)
inModule = is.finite(match(moduleColors, modules));
modProbes = probes[inModule];
#modGenes = annot$gene_symbol[match(modProbes, annot$substanceBXH)];
# Select the corresponding Topological Overlap
modTOM = TOM[inModule, inModule];
dimnames(modTOM) = list(modProbes, modProbes)
# Export the network into edge and node list files Cytoscape can read
cyt = exportNetworkToCytoscape(modTOM,
                               edgeFile = paste("AS-green-FPKM-Step-by-step-CytoscapeInput-edges-", paste(modules, collapse="-"), ".txt", sep=""),
                               nodeFile = paste("AS-green-FPKM-Step-by-step-CytoscapeInput-nodes-", paste(modules, collapse="-"), ".txt", sep=""),
                               weighted = TRUE,
                               threshold = 0.02,
                               nodeNames = modProbes,
                               #altNodeNames = modGenes,
                               nodeAttr = moduleColors[inModule]);
#=====================================================================================
#  分析网络可视化,用heatmap可视化权重网络,heatmap每一行或列对应一个基因,颜色越深表示有较高的邻近
#=====================================================================================
options(stringsAsFactors = FALSE);
lnames = load(file = "AS-green-FPKM-01-dataInput.RData");
lnames
lnames = load(file = "AS-green-FPKM-02-networkConstruction-stepByStep.RData");
lnames
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
#1. 可视化全部基因网络
# Calculate topological overlap anew: this could be done more efficiently by saving the TOM
# calculated during module detection, but let us do it again here.
dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6);
# Transform dissTOM with a power to make moderately strong connections more visible in the heatmap
plotTOM = dissTOM^7;
# Set diagonal to NA for a nicer plot
diag(plotTOM) = NA;
# Call the plot function
#sizeGrWindow(9,9)
TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
#随便选取1000个基因来可视化
nSelect = 1000
# For reproducibility, we set the random seed
set.seed(10);
select = sample(nGenes, size = nSelect);
selectTOM = dissTOM[select, select];
# There's no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select];
# Open a graphical window
#sizeGrWindow(9,9)
# Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing
# the color palette; setting the diagonal to NA also improves the clarity of the plot
plotDiss = selectTOM^7;
diag(plotDiss) = NA;
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")


根据基因间表达量进行模块聚类(前边得到)和所得到的各模块间的相关性

MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
MET = orderMEs(MEs)
sizeGrWindow(7, 6) 
plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2), plotDendrograms = FALSE, xLabelsAngle = 90)


脚本来自:
http://tiramisutes.github.io/2016/09/14/WGCNA.html

https://bioconductor.org/packages/devel/bioc/vignettes/CVE/inst/doc/WGCNA_from_TCGA_RNAseq.html

http://www.bio-info-trainee.com/2535.html

WGCNA原理及应用

WGCNA介绍:

WGCNA(weighted gene co-expression network analysis,权重基因共表达网络分析)是一种分析多个样本基因表达模式的分析方法,可将表达模式相似的基因进行聚类,并分析模块与特定性状或表型之间的关联关系,因此在疾病以及其他性状与基因关联分析等方面的研究中被广泛应用。

WGCNA算法是构建基因共表达网络的常用算法(详解:http://www.jianshu.com/p/94b11358b3f3)。WGCNA算法首先假定基因网络服从无尺度分布,并定义基因共表达相关矩阵、基因网络形成的邻接函数,然后计算不同节点的相异系数,并据此构建分层聚类树(hierarchical clustering tree),该聚类树的不同分支代表不同的基因模块(module),模块内基因共表达程度高,而分属不同模块的基因共表达程度低。最后,探索模块与特定表型或疾病的关联关系,最终达到鉴定疾病治疗的靶点基因、基因网络的目的。在该方法中module被定义为一组具有类似表达谱的基因,如果某些基因在一个生理过程或不同组织中总是具有相类似的表达变化,那么我们有理由认为这些基因在功能上是相关的,可以把他们定义为一个模块(module)。这似乎有点类似于进行聚类分析所得到结果,但不同的是,WGCNA的聚类准则具有生物学意义,而非常规的聚类方法(如利用数据间的几何距离),因此该方法所得出的结果具有更高的可信度。当基因module被定义出来后,我们可以利用这些结果做很多进一步的工作,如关联性状,代谢通路建模,建立基因互作网络等。

 

WGCNA的用处:

 

  • 这类处于调控网络中心的基因称为核心基因(hub gene),这类基因通常是转录因子等关键的调控因子,是值得我们优先深入分析和挖掘的对象。
  • 在网络中,被调控线连接的基因,其表达模式是相似的。那么它们潜在有相似的功能。所以,在这个网络中,如果线条一端的基因功能是已知的,那么就可以预测线条另一端的功能未知的基因也有相似的功能。

 

下面的问答来自基迪奥,也能加深对WGCNA的理解

问1、调控网络和共表达网络有什么区别?
答:调控网络是个更广泛的概念,而共表达网络是调控网络的一种。
理论上我们可以利用各类信息构建调控网络(表达相关性,序列靶向关系、蛋白互作关系),另外调控网络构建的信息既可以来源真实的实验验证的关系,也可以来源生物信息的预测。而共表达网络特指利用基因间的表达相关性预测基因间调控关系的方法,而WGCNA又是共表达网络分析中最有效的方法之一。

问2、WGCNA分析适合的生物物种范围有规定么?
答:没有限制。对于任何物种中心法则都是存在的,调控关系对于任何物种都是存在的,所以WGCNA没有物种限定。

问3、同一物种,不同来源的转录组数据(比如不同文章/资料来源的),可以放在一起做WGCNA分析吗?
答:只要样本间有相似的生物学意义,是可以合并在一起做分析的。但要注意,不同批次之间的样本是有批次效应的,所以可能会带来一些误差,但是是可以放在一起分析的。

问4、相同材料不同处理之间,可以放在一起做WGCNA分析吗?比如重金属和盐碱处理。
答:可以的。这也正式WGCNA强大的地方,其可以将不同处理的样本,合并在一起做分析。其他方法则不一定有这么强大的能力,比如做基因表达趋势分析时,如果样本涉及到多个处理不同时期的时候,就不好合并分析(或合并后难以解读)。但WGCNA的方法关注的是调控关系,所以不管是多少个处理组,都可以很好的整合在一起做分析。

问5、不同批次的数据能放一起做WGCNA吗?
答:可以的。虽然有批次的干扰,但是干扰对WGCNA网络没有太大影响。因为WGCNA不是做差异分析,而是基因的共表达。因为批次效应理论上不影响相关性。

问6、不同类型的材料,比如亲本和F1,适合放一起进行WGCNA么?
答:如果是一个作图群体,当然亲本与F1是可以放在一起分析的,因为你只关心基因的表达模式,所以把亲本加进来是没有问题的。

问7、没有生物学重复,共3组,每组5个时间点能够做吗?
答:理论上有15个样本,是可以做WGCNA分析的。并且,分析出来的结果对你的研究应该是非常有用的。至少他会比趋势分析更有意义,更加准确。

问8、一般说WGCNA的样品不少于15个,15个样品考虑重复吗?不同倍性的材料呢?
答:15个样本这个是包含了生物学重复,比如5个时间点3个重复;在RNA-seq里面建议不要用不同倍性材料加进来。除非是有参考的多倍体,如果是无参的多倍体,不同倍性之间差异太大,会让调控网络不准确。所以用单一倍性的材料做调控网络会更加准确。

问9、可以将RNA-seq数据与蛋白组数据,甲基化数据放一起做WGCNA分析?
答:不能与蛋白数据一起分析。因为WGCNA是基于相关系数的算法。所以最好一起分析的数据变异度是类似的,RNAseq变异非常大,而蛋白的数据变异很小,两者的变化不在一个数量级上面。所以两种数据放在一起分析不合理。
但RNA数据可以尝试跟甲基化数据一起分析。当然我们也建议RNA数据与代谢组数据一起分析,因为代谢组的数据变异也非常大。

问10、表达量和表达的基因数目差异太大的样品可以一起分析吗?比如样品A有2k个gene表达 而样品B有2w个gene表达了 AB可以一起分析吗?
答:做WGCNA分析的时候,不能脱离生物学意义,既然要分析调控网络,那么应该分析有相似生物学意义的一组基因,比如说拿相似组织来一起做分析,比如不应该拿大脑的样本与脚趾的样本合并在一起做分析,因为很显然,这两个组织没有关联。如果两个样本之间是有相关联的生物学意义,哪怕表达的基因数不一样,或表达模式差异很大,那依然可以放在一起分析;但如果样本之间完全没有生物学意义,那么分析就没有意义。

问11、实验设计是case3个时间点(各点都有三个重复),control同样的3个时间点(每点三个重复),WGCNA怎么做?3个时间点和case-control两个因素能同时考虑进来分析吗?
答:可以的。做WGCNA是更加合理的,因为有两个梯度的样本,如果只是做差异分析的话,逻辑可能非常复杂,做WGCNA分析是对样本特性更好的解析,可以直观看到基因在六个处理组里面是怎样表达的。

问12、可以拿混合样本分析吗?比如一个病原细菌跟人类细胞的基因,能说明细菌跟人类细胞基因有调控关系吗?
答:可以。前提是病原菌有足够的数据并定量准确,并且这个分析是非常有意义的,最后可以说明这些病原菌可以调控哪些宿主基因。

问13、但是病原宿主混合分析的话,宿主蛋白不能分泌到宿主体内岂不是WGCNA生物学上也没有意义吗?
答:依然有意义。即使病原的基因没有分泌到宿主里面,但是病原的蛋白是会影响宿主基因的调控的,比如某个细菌感染某个植物,虽然细菌的蛋白不能直接分泌到植物体内,但会影响植物蛋白的分泌。混在一起分析依然是有意义,可以看到植物里面到底哪个基因对细菌蛋白产生应答作用。

问14、芯片数据两分类,每组20个样本,能否每组单独做WGCNA?
答:可以。WGCNA还有一种重要功能是做两个网络的比较,比如病人20个样本做一个调控网络,健康人做一个调控网络,然后两个网络做比较。

问15、WGCNA可以用来分析lncRNA对下游基因的调控分析吗?
答:可以。WGCNA网络有利于预测lncRNA的潜在功能。

问16、构建网络是用所有表达基因还是差异基因?
答:这个是具体问题具体分析。如果使用所有的基因分析,会导致运算量非常大。而也不是所有的基因在这个实验中都有生物学意义,所以我们会提前做一些过滤。
但用于分析的基因不一定是差异表达基因,有时可以用差异表达基因做一个并集,或通过计算变异系数将变异系数低的基因以及低表达的基因去除。但注意,如果你有关心的特定目标基因的话,应该尽量给予保留。

问17、关注某一个pathway上的基因以及调控因子之间的相关性,构建WGCNA网络的时候属于这个pathway的基因数量太少会不会影响结果呢?
答:这不是问题。在一个调控网络里面,样本的某个pathway上,并不是所有基因参与调控(或存在差异性),所以在做WGCNA分析的时候,会做一些过滤,将有变化的基因挑出来再做分析。即分析的是某个pathway上有变化的基因,不需要分析pathway上所有的基因,只需要分析那些变化的基因就够了。

问18、前期筛选的时候,要选出在所有样本中变异系数比较大的基因呢?还是直接用差异表达的基因取并集?用基因还是转录本,哪个好呢?
答:两则都可以,我推荐使用变异系数,选择那些变异较大的基因,来做下面的分析。然后建议用基因不要用转录本,因为转录本的定量是不准确的。

问19、变异系数一般取多大?

答:具体问题具体分析。例如,没有特定目标的时候,可以先计算变异系数,将变异系数的百分之前50来做分析,把变异系数偏低的后面一半过滤掉。

问20、输入数据用FPKM合适吗?
答:可以。

问21、RNA seq数据是RSEM值怎么办?
答:RSEM值原始输出结果为reads数,如果是RSEM值建议做一个RPKM校正再做分析。

问22、除了RPKM值以外,做WGANA是否还需要其他数据?TCGA数据可否来做WGCNA分析?
答:在做WGCNA分析必须要用表达量数据,但TCGA的数据某些层级没有表达量数据,没有表达量数据自然就无法做WGCNA分析。

问23、请问输入的基因样本的矩阵的时候,要不要对数据标准化?
答:做WGCNA分析的时候,不需要对数据进行标准化,输入RPKM值就足以做这个分析。虽然一些文章会做log2处理,但我认为取了LOG2后,会让一些表达关系没有那么丰富。

问24、每个样本有3个生物学重复,不需要对三个重复的表达量求平均值代表该样本吗?
答:注意,做WGCNA的时候每个样本是独立的,三个生物学重复样本是全部导入做分析,不是取均值再做分析,每个样本都是独立的。

问25、如果3个生物学重复,做WGCNA的时候是取三个值,还是用cuffdiff处理后取一个值?
答:如果是生物学重复样本进行调控网络分析,每个样本独立使用,而不是取均值。

问26、请问将样本信息同模块特征值进行相关性分析的时候,样本信息是怎么处理的呢?比如不同取样点、不同性别什么的,这不是数量性状信息的,这种情况应该怎么处理呢?
答:样本的任何信息都可以做模块相关性分析。比如相关时间点,可以按照先后量化为12134567。又如不同性别,男与女,可以定义为1,-1。任何性状量化为数字后,都可以进行相关性分析。

问27、怎么将模块与性状对应起来呢有些性状不好量化,如果直接将模块与分组对应,如何实现, 不需要量化指标么?
答:首先需要将性状量化,如果无法将性状量化,那么就无法分析。至于分组信息,也可以量化为类似00001111000(1代表一种组别,2代表另一组组别),实现分组信息的数字化。

问28、基因数量为3w左右时,modules数量为多少结果较为理想?怎么评价聚类效果的好坏?
答:modules数量没有标准,modules数量无法评估模块分的好坏,分组是否合理应该看树的树形图,比如树的分支很清晰就说明模块式清晰的。modules数量数由生物性状决定的。比如样本表达信息很丰富的时候,modules数量会很多;如果样本的基因表达相对单一,modules数量就会比较少。

问29、我运行例子的时候,得出来基因之间的direction全是undirected,这和前面的几种关系有什么区别?
答:WGCNA是一个undirected的方法,它的网络是无方向的,有相关关系但是无方向。

问30、如果做有向网络的构建,您推荐那些方法?
答:很多方法,例如贝叶斯的方法。

问31、非模式物种可以得出基因之间的相互关系类型么?得出的结果也是undirected么?
答:WGCNA是基于表达两处理的,所以即使是非模式生物,当然也可以他们之间关系,并且关系也是一个无向网络。

问32、选择几个表型数据进行结合分析比较好
答:越多越好,看实验设计。

问33、感染小鼠,5个时间点,3个重复,找不到合适的表型怎么办?
答:如果找不到合适表型,可以找某个时间点应答的基因,本身基因的表达趋势已经有某种生物学意义的。没有找到合适表型,也可以看变化趋势。不一定要做表型的相关分析,其他分析也是很有趣的。例如,可以对模块功能的富集分析,其实都是可以帮助你找到特定模块的。所以不用纠结于做某个表型的关联分析。

问34、weight就是tom值吗?
答:是的。

问35、剪模块是怎么做的?是根据TOM划分吗?需要自己设定,还是R自动的?
答:剪模块是R中自动完成的,不需要划分,但合并的时候你可以设定一个指标,比如差异度是0.25。

问36、看WGCNA说明是用相异矩阵D(D=1-TOM)去做聚类,然后动态剪切?
答:用TOM值来构建矩阵,TOM值就是两个样本的相似度,1-TOM值就是两个样本的差异度,相似度与差异度可以理解为一个东西,并不矛盾。

问37、模块特征值和样本性状相关分析的具体方法是?
答:R包用的是计算相关系数的方法。

问38、WGCNA里面一般会提到hubgene,如何确定hubgene?
答:在WGCNA分析里面,每个基因都会计算连通性,连通性高的就是hubgene。

问39、在R中安装“”WGCNA“”说不适合R3.3.1,那适合哪个版本?
答:WGCNA应该是所有版本都适合,如果版本没有可以考虑降低R软件的版本,这个对分析没有影响。因为不同R版本是一样的。

问40、用STEM分析的时候拟合多少个模型合适?
答:建议不要超过20个。模块太多不好分析。

参考网站:

http://tiramisutes.github.io/2016/09/14/WGCNA.html

http://www.jianshu.com/p/94b11358b3f3

http://www.omicshare.com/class/home/index/classdetail?id=20