pipeline: QC >>> alignment >>> differential expression analysis >>> enrichment analysis

primary softs or packages involved in the pipeline：

1. QC (linux)

• fastqc: quality assessment
• cutadapt: remove Illumina Universal Adaptor or overrepresentative sequences
• fastx_trimmer: trim the low-quality bases (e.g the first 10bp)
• multiqc: intergrate results of fastqc

2. alignment (linux)

• STAR: faster, ~90 minutes for each sample
• bwa

3. Differential expression analysis (R packages)

• HTseq: calculate raw read counts for each transcript or gene
• DESeq2
• limma

4. enrichment analysis

• clusterProfiler

5. plot

• ggplot2
• Cairo
• pheatmap
• RColorBrewer

main codes
input: raw counts
output: differential expression analysis and enrichment analysis based on DEGs
输出差异分析结果列表、火山图、基于差异表达基因的富集分析结果(GO/KEGG/ReactomePA)列表,适合粗略地先看下结果。

Volcano plot

#载入R包
require("DESeq2")
require("limma")
require("pasilla")
require("vsn")
require("ggplot2")
require("Cairo")
require("pheatmap")
require("RColorBrewer")
require(DOSE)
require(org.Mm.eg.db)
require(clusterProfiler)
require(stringr)
require(ReactomePA)
source('deg_function.r') #cal DEGs
source('enrich_function.r') #enrichment analysis

#读入raw count
d<-read.table("combined.count",row.names=1,header = T)
head(d,2)
#                    A1  A2  A3  B1  B2  B3  C1  C2  C3  D1  D2  D3  
#ENSMUSG00000000001 615 645 531 767 684 686 751 778 929 767 680 785
#ENSMUSG00000000003 0   0   0   0   0   0   0   0   0   0   0   0

tail(d)
#                    A1  A2  A3  B1  B2  B3  C1  C2  C3  D1  D2  D3  
#ENSMUSG00000116528 0   0   0   0   0   0   0   0   0   0   0   0
#__no_feature   2509484 2534451 1672031 2405219 2323971 2398776 2361692 2600649 2961624 2412725 2331754 2562243
#__ambiguous    1198166 1210490 1057190 1315088 1268587 1281891 1225386 1328553 1443758 1292656 1255907 1234649
#__too_low_aQual    0   0   0   0   0   0   0   0   0   0   0   0
#__not_aligned  0   0   0   0   0   0   0   0   0   0   0   0
#__alignment_not_unique 796927  831108  756030  936088  901476  884335  862996  827094  976800  856331  789099  848550

dat<-d[-c(53801:53806),] #删除最后几行信息

#ID convert:ENSEMBL to SYMBOL
ee=bitr(rownames(dat),fromType = "ENSEMBL",toType = "SYMBOL", OrgDb="org.Mm.eg.db")
write.csv(ee,file="ensembl_symbol.csv",quote=F)

#set groups
a=c(1:3);b=c(4:6);c=c(7:9);d=c(10:12);
caselist=list(a,b,c,d) 
namelist=c("A","B","C","D") 

#all samples included: is there any outlier?
a=c(rep("case",9),rep("ctrl",3)) #前9列为case，最后3列为ctrl
name="all"
dd=dat
deg_function(a,dd,name)

#comparison among groups: A vs. B/C/D; B vs C/D; C vs D (the first group serves as controls)
for (i in 1:3) {for(j in (i+1):4){
    a=c(rep("ctrl",3),rep("case",3))
    name=paste(namelist[i],namelist[j],sep="_")
    dd=dat[,c(caselist[[i]],caselist[[j]])]
    deg_function(a,dd,name)    
    } 
}

deg_function.r: 这里仅取p adj<0.05 & abs(log2FC)>1的基因做富集分析，如需调整参数，可修改up和down

deg_function <- function(a,dd,name){
  condition=factor(a)
  coldata<-data.frame(condition)
  rownames(coldata)<-colnames(dd)
  ds<-DESeqDataSetFromMatrix(countData=dd,colData=coldata,design = ~ condition)
  #pre-filtering
  keep <- rowSums(counts(ds))>=10
  ds <- ds[keep,]
  nrow(ds) 
  ds$condition <- relevel(ds$condition,ref="ctrl") #指定对照组
  #normalization后计算样本间的距离:sampe vs. sample
  rld <- rlog(ds, blind=FALSE)
  sampleDists <- dist(t(assay(rld)))
  sampleDistMatrix <- as.matrix(sampleDists)
  rownames(sampleDistMatrix) <- colnames(rld)
  colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
  #样本间相似性图
  CairoPNG(paste(name,'distance.png',sep=""))
  pheatmap(sampleDistMatrix,
           clustering_distance_rows=sampleDists,
           clustering_distance_cols=sampleDists,
           col=colors)
  dev.off()
  
  #输出归一化后的count数据
  dds <- estimateSizeFactors(ds)
  norcounts <- counts(dds, normalized=T)
  write.csv(norcounts,file=paste(name,"norm_counts.csv",sep=""),quote=F)
  #DEG
  des <- DESeq(ds)
  res <- results(des)
  resOrdered <- res[order(res$padj),]
  resOrdered_out=cbind(resOrdered@rownames,as.data.frame(resOrdered))
  colnames(resOrdered_out)[1]=colnames(ee)[1]
  out=merge(resOrdered_out,ee,by="ENSEMBL",all.x=TRUE)  
  write.csv(out,file=paste(name,"_DEG.csv",sep=""),quote = F)  
  #出火山图
  CairoPNG(paste(name,'Volcano.png',sep=""))
  with(resOrdered, plot(resOrdered$log2FoldChange, -log10(resOrdered$padj), pch=1, main="Volcano plot",xlab="log2FC",ylab="-log10(padj)"))
  with(subset(resOrdered, padj<0.05 ), points(log2FoldChange, -log10(padj), pch=1, col="red"))
  with(subset(resOrdered, abs(log2FoldChange)>1.5), points(log2FoldChange, -log10(padj), pch=1, col="orange"))
  with(subset(resOrdered, padj<0.05 & abs(log2FoldChange)>1.5), points(log2FoldChange, -log10(padj), pch=1, col="green"))
  dev.off()
  
  #enrichment analysis
  up=subset(resOrdered,padj<0.05 & log2FoldChange>1)
  x1=rownames(up);y1=paste(name,"up",sep="_")
  if(length(x1)>5){enrich_function(x1,y1)}  #DEG数目过少，没有mapped ENTRE ID的话，会报错中断
  down=subset(resOrdered,padj<0.05 & log2FoldChange < -1)
  x2=rownames(down);y2=paste(name,"down",sep="_")
  if(length(x2)>5){enrich_function(x2,y2)}
}

enrich_function.r: 现设置为mouse，换物种需修改参数 OrgDb

enrich_function<-function(x,y){
  
  eg = bitr(x, fromType="ENSEMBL", toType="ENTREZID", OrgDb="org.Mm.eg.db")
  gene = eg$ENTREZID
  
  BP <- enrichGO(gene, "org.Mm.eg.db", keyType = "ENTREZID",ont = "BP",pvalueCutoff  = 0.05,pAdjustMethod = "BH",qvalueCutoff  = 0.1, readable=T)
  MF <- enrichGO(gene, "org.Mm.eg.db", keyType = "ENTREZID",ont = "MF",pvalueCutoff  = 0.05,pAdjustMethod = "BH",qvalueCutoff  = 0.1, readable=T)
  CC <- enrichGO(gene, "org.Mm.eg.db", keyType = "ENTREZID",ont = "CC",pvalueCutoff  = 0.05,pAdjustMethod = "BH",qvalueCutoff  = 0.1, readable=T)
  BP_simp <- simplify(BP, cutoff=0.1,by="p.adjust",select_fun=min)
  MF_simp <- simplify(MF, cutoff=0.1,by="p.adjust",select_fun=min)
  CC_simp <- simplify(CC, cutoff=0.1,by="p.adjust",select_fun=min)
  bp=cbind(rep("BP",times=nrow(as.data.frame(BP_simp@result))),as.data.frame(BP_simp@result))
  mf=cbind(rep("MF",times=nrow(as.data.frame(MF_simp@result))),as.data.frame(MF_simp@result))
  cc=cbind(rep("CC",times=nrow(as.data.frame(CC_simp@result))),as.data.frame(CC_simp@result))
  colnames(cc)=colnames(bp)
  colnames(mf)=colnames(bp)
  write.csv(rbind(bp,mf,cc),file=paste(y,"GO.csv",sep="_"),quote=F)
  
  kk <- enrichKEGG(gene = gene,
                   organism ='mmu',
                   pvalueCutoff = 0.05,
                   qvalueCutoff = 0.1,
                   minGSSize = 1,
                   #readable = TRUE ,
                   use_internal_data =FALSE)
  write.csv(as.data.frame(kk@result), file=paste(y,"kegg.csv",sep="_"),quote=FALSE)
  
  react <- enrichPathway(gene=gene,pvalueCutoff=0.05, readable=T,organism = "mouse")
  write.csv(react,file=paste(y,"ReactomePA.csv",sep="_"),quote=F)
}