• 前面RNA-seq分析：从软件安装到富集分析部分已经把转录组全部流程走完了一遍，这次利用RNA-seq(2)-2:下载数据中下载的肝癌数据进行分析，不在赘述细节，所以有看细节的还是请去这里。
• 这部分主要是代码和部分结果图，但进行了部分修正，比如KEGG 可视化部分，用了前clusterprofiler部分的结果，所以这部分不包括Gage包。
• 所以从质量比对开始进行

-----------------------------------------------------------------

3 sra到fastq格式转换并进行质量控制

3.1 数据解压：用samtools中的fastq-dump将sra格式转为fastq格式

• 备注：用时5个小时

for ((i=2;i<=5;i++));do fastq-dump --gzip --split-3 -A SRR31621$i.sra -O .;done

3.2 用fastqc进行质量控制

#将所有的数据进行质控，得到zip的压缩文件和html文件
fastqc -o .  *.fastq.gz

3.3 3 质控结果解读(为保持紧凑，只放一张图）

SRR316214.png

4 下载参考基因组及基因注释

RNA-seq(4):下载参考基因组及基因注释部分已经下载

5 序列比对：Hisat2

5.1 开始比对：用hisat2，得到SAM文件(5个小时)

• 我的fastq文件在/mnt/f/rna_seq/data
• 我的reference在/mnt/f/rna_seq/data/reference
• 我的index在/mnt/f/rna_seq/data/reference/index/hg19
• 比对后得到的bam文件会存放在/mnt/f/rna_seq/aligned

 source ~/miniconda3/bin/activate
(base) kelly@DESKTOP-MRA1M1F:/mnt/f/rna_seq/data$ cd ..
(base) kelly@DESKTOP-MRA1M1F:/mnt/f/rna_seq$
for ((i=2;i<=5;i++));do hisat2 -t -x /mnt/f/rna_seq/data/reference/index/hg19/genome -1 /mnt/f/rna_seq/data/SRR31621${i}.sra_1.fastq.gz -2 /mnt/f/rna_seq/data/SRR31621${i}.sra_2.fastq.gz -S SRR31621${i}.sam ;done

(base) kelly@DESKTOP-MRA1M1F:/mnt/f/rna_seq$ for ((i=2;i<=5;i++));do hisat2 -t -x /mnt/f/rna_seq/data/reference/index/hg19/genome -1 /mnt/f/rna_seq/data/SRR31621${i}.sra_1.fastq.gz -2 /mnt/f/rna_seq/data/SRR31621${i}.sra_2.fastq.gz -S SRR31621${i}.sam ;done
Time loading forward index: 00:00:05
Time loading reference: 00:00:01
Multiseed full-index search: 02:11:39
101339131 reads; of these:
  101339131 (100.00%) were paired; of these:
    7840870 (7.74%) aligned concordantly 0 times
    87196406 (86.04%) aligned concordantly exactly 1 time
    6301855 (6.22%) aligned concordantly >1 times
    ----
    7840870 pairs aligned concordantly 0 times; of these:
      333400 (4.25%) aligned discordantly 1 time
    ----
    7507470 pairs aligned 0 times concordantly or discordantly; of these:
      15014940 mates make up the pairs; of these:
        8711709 (58.02%) aligned 0 times
        5441734 (36.24%) aligned exactly 1 time
        861497 (5.74%) aligned >1 times
95.70% overall alignment rate

5.2 SAM文件转换为BAM文件,并对bam文件进行sort,最后建立索引：SAMtools

Sam 文件转换为bam格式（用时1.5hrs）

(base) kelly@DESKTOP-MRA1M1F:/mnt/f/rna_seq/aligned$ for ((i=2;i<=5;i++));do samtools view -S /mnt/f/rna_seq/SRR31621${i}.sam -b > SRR31621${i}.bam;done

对bam文件进行排序，默认染色体位置（用时1.5hrs）

for ((i=2;i<=5;i++));do samtools sort SRR31621${i}.bam
-o SRR31621${i}_sorted.bam;done

建立索引（用时40mins）

for ((i=2;i<=5;i++));do samtools index SRR31621${i}_sorted.bam;done

6 reads计数，合并矩阵并进行注释

6.1 bam文件按reads name排序（用时1h）

for ((i=2;i<=5;i++));do samtools sort -n SRR31621${i}.b
am -o SRR31621${i}_nsorted.bam;done

6.2 2 reads计数，得到表达矩阵htseq-count（用时很久很久，久到不忍写，8hrs）

注释文件如果已经解压，则不再需要重复

cd ../data/matrix
gunzip /mnt/f/rna_seq/data/reference/annotation/hg19/gencode.v19.annotation.gt.gz && rm -rf gencode.v19.annotation.gtf.gz
for ((i=2;i<=5;i++));do htseq-count -r name -f bam /mnt/f/rna_seq/aligned/SRR31621${i}_nsorted.bam /mnt/f/rna_seq/data/reference/annotation/hg19/gencode.v19.an
notation.gtf > SRR31621${i}.count; done

ls -al *.count
-rwxrwxrwx 1 root root 1197426 Aug  7 16:25 SRR316212.count
-rwxrwxrwx 1 root root 1186189 Aug  7 17:16 SRR316213.count
-rwxrwxrwx 1 root root 1200305 Aug  7 22:51 SRR316214.count
-rwxrwxrwx 1 root root 1187596 Aug  7 23:32 SRR316215.count

6.3 3 合并表达矩阵并进行基因名注释（R中进行）

Tumor:SRR316214,SRR316215
Adjacent Normal Liver:SRR316212,SRR316213

setwd("F:/rna_seq/data/matrix")
options(stringsAsFactors = FALSE)
control1<-read.table("SRR316212.count",sep = "\t",col.names = c("gene_id","control1"))
control2<-read.table("SRR316213.count",sep = "\t",col.names = c("gene_id","control2"))
treat1<-read.table("SRR316214.count",sep = "\t",col.names = c("gene_id","treat1"))
treat2<-read.table("SRR316215.count",sep = "\t",col.names = c("gene_id","treat2"))
raw_count <- merge(merge(control1, control2, by="gene_id"), merge(treat1, treat2, by="gene_id"))
head(raw_count)
raw_count_filt <- raw_count[-1:-5,]
ENSEMBL <- gsub("\\.\\d*", "", raw_count_filt$gene_id) 
row.names(raw_count_filt) <- ENSEMBL
head(raw_count_filt)

> head(raw_count_filt)
                           gene_id control1 control2 treat2.x treat2.y
ENSG00000000003 ENSG00000000003.10     4004      781      756      756
ENSG00000000005  ENSG00000000005.5        1        0        0        0
ENSG00000000419  ENSG00000000419.8      776      140      165      165
ENSG00000000457  ENSG00000000457.9      624      144      240      240
ENSG00000000460 ENSG00000000460.12      260       52      105      105
ENSG00000000938  ENSG00000000938.8      161       59       16       16

7 DEseq2筛选差异表达基因并注释

这次我换了Annotation包进行注释

7.1 载入数据（countData和colData）

# 这一步很关键，要明白condition这里是因子，不是样本名称；小鼠数据有对照组和处理组，各两个重复
condition <- factor(c(rep("control",2),rep("treat",2)), levels = c("control","treat"))
condition
colData <- data.frame(row.names=colnames(mycounts), condition)

> colData
         condition
control1   control
control2   control
treat1       treat
treat2       treat

7.2 构建dds对象,开始DESeq流程

dds <- DESeqDataSetFromMatrix(countData=mycounts, 
                              colData=colData, 
                              design= ~ condition)
dds = DESeq(dds)

7.3 总体结果查看

接下来，我们要查看treat versus control的总体结果，并根据p-value进行重新排序。利用summary命令统计显示一共多少个genes上调和下调（FDR0.1）

res = results(dds, contrast=c("condition", "control", "treat"))
res = res[order(res$pvalue),]
head(res)
summary(res)
write.csv(res,file="All_results.csv")
table(res$padj<0.01)

> summary(res)
out of 30981 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)     : 3928, 13% 
LFC < 0 (down)   : 3283, 11% 
outliers [1]     : 0, 0% 
low counts [2]   : 10317, 33% 
(mean count < 2)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
> table(res$padj<0.01)
FALSE  TRUE 
16192  4472

可见，上调基因和下调的数目都很多，padj<0.01的有4472，即使padj<0.001的也有3089个。

7.4 提取差异表达genes（DEGs）

先看下padj（p值经过多重校验校正后的值）小于0.01，表达倍数取以2为对数后大于1或者小于-1的差异表达基因。代码如下

diff_gene_deseq2 <-subset(res, padj < 0.01 & abs(log2FoldChange) > 1)
dim(diff_gene_deseq2)
head(diff_gene_deseq2)
write.csv(diff_gene_deseq2,file= "HCC_DEG_0.05_2.csv")

> dim(diff_gene_deseq2)
[1] 3780    6
> head(diff_gene_deseq2)
log2 fold change (MLE): condition control vs treat 
Wald test p-value: condition control vs treat 
DataFrame with 6 rows and 6 columns
                 baseMean log2FoldChange     lfcSE      stat        pvalue          padj
                <numeric>      <numeric> <numeric> <numeric>     <numeric>     <numeric>
ENSG00000130649 78059.342       7.190959 0.2110024  34.07998 1.460506e-254 3.017990e-250
ENSG00000268925  5666.055       7.534595 0.2573299  29.27991 1.869559e-188 1.931628e-184
ENSG00000105697  8791.483       6.611356 0.2338215  28.27523 6.970959e-176 4.801597e-172
ENSG00000167244  8988.099       6.767261 0.2582403  26.20528 2.313301e-151 1.195051e-147
ENSG00000162366  4037.960      -7.076854 0.2704122 -26.17062 5.742578e-151 2.373293e-147
ENSG00000169715  3803.597       6.238596 0.2403373  25.95767 1.489813e-148 5.130914e-145

把padj<0.001,|log2FC|>1有2938个

8 探索分析结果:Data visulization

#MA plot
plotMA(res,ylim=c(-2,2))
topGene <- rownames(res)[which.min(res$padj)]
with(res[topGene, ], {
  points(baseMean, log2FoldChange, col="dodgerblue", cex=6, lwd=2)
  text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
})

#shirnked
res_order<-res[order(row.names(res)),]
res = res_order
res.shrink <- lfcShrink(dds, contrast = c("condition","treat","control"), res=res)
plotMA(res.shrink, ylim = c(-5,5))
topGene <- rownames(res)[which.min(res$padj)]
with(res[topGene, ], {
  points(baseMean, log2FoldChange, col="dodgerblue", cex=2, lwd=2)
  text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
})

#identify
idx <- identify(res$baseMean, res$log2FoldChange)
rownames(res)[idx]
#plotcounts
plotCounts(dds, gene=which.min(res$padj), intgroup="condition", returnData=TRUE)
plotCounts(dds, gene="ENSG00000130649", intgroup="condition", returnData=FALSE)
#boxplot
plotCounts(dds, gene="ENSG00000130649", intgroup="condition", returnData=TRUE) %>% 
  ggplot(aes(condition, count)) + geom_boxplot(aes(fill=condition)) + scale_y_log10() + ggtitle("CYP2E1")
#pointplot
d <- plotCounts(dds, gene="ENSG00000130649", intgroup="condition", returnData=TRUE)
ggplot(d, aes(x=condition, y=count)) + 
  geom_point(aes(color= condition),size= 4, position=position_jitter(w=0.5,h=0)) + 
  scale_y_log10(breaks=c(25,100,400))+ ggtitle("CYP2E1")

##3最小padj
d <- plotCounts(dds, gene=which.min(res$padj), intgroup="condition", 
                returnData=TRUE)
ggplot(d, aes(x=condition, y=count)) + 
  geom_point(position=position_jitter(w=0.1,h=0)) + 
  scale_y_log10(breaks=c(25,100,400))

#PCA
vsdata <- vst(dds, blind=FALSE)
plotPCA(vsdata, intgroup="condition")
#beatifule pca plot
pcaData <- plotPCA(vsdata, intgroup=c("condition"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(PC1, PC2, color=condition, shape=condition)) +
  geom_point(size=3) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  coord_fixed()

library("pheatmap")
select<-order(rowMeans(counts(dds, normalized = TRUE)),
              decreasing = TRUE)[1:20]
df <- as.data.frame(colData(dds)[,c("condition","sizeFactor")])
# this gives log2(n + 1)
ntd <- normTransform(dds)
pheatmap(assay(ntd)[select,], cluster_rows=FALSE, show_rownames=FALSE,
         cluster_cols=FALSE, annotation_col=df)
#vsd
vsd <- vst(dds, blind=FALSE)
pheatmap(assay(vsd)[select,], cluster_rows=FALSE, show_rownames=FALSE,
         cluster_cols=FALSE, annotation_col=df)
##sample to sample heatmap
sampleDists <- dist(t(assay(vsd)))
library("RColorBrewer")
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- paste(vsd$condition, vsd$type, sep="-")
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
         clustering_distance_rows=sampleDists,
         clustering_distance_cols=sampleDists,
         col=colors)

这部分主要结果部分的一些可视化，具体参考详细的部分即可，只放一张图

Rplot.jpeg

9 DEGs的富集分析(功能注释）ClusterProfiler包

##enrichment analysis using clusterprofiler package created by yuguangchuang
library(clusterProfiler)
library(DOSE)
library(stringr)
library(org.Hs.eg.db)
#get the ENTREZID for the next analysis
sig.gene= diff_gene_deseq2
head(sig.gene)
gene<-rownames(sig.gene)
head(gene)
gene.df<-bitr(gene, fromType = "ENSEMBL", 
              toType = c("SYMBOL","ENTREZID"),
              OrgDb = org.Hs.eg.db)

head(gene.df)

> head(gene.df)
          ENSEMBL   SYMBOL ENTREZID
1 ENSG00000130649   CYP2E1     1571
3 ENSG00000105697     HAMP    57817
4 ENSG00000167244     IGF2     3481
5 ENSG00000162366 PDZK1IP1    10158
6 ENSG00000169715     MT1E     4493
7 ENSG00000074276    CDHR2    54825

GO enrichment

#Go enrichment
ego_cc<-enrichGO(gene       = gene.df$ENSEMBL,
                 OrgDb      = org.Hs.eg.db,
                 keyType    = 'ENSEMBL',
                 ont        = "CC",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.01,
                 qvalueCutoff = 0.05)
ego_bp<-enrichGO(gene       = gene.df$ENSEMBL,
                 OrgDb      = org.Hs.eg.db,
                 keyType    = 'ENSEMBL',
                 ont        = "BP",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.01,
                 qvalueCutoff = 0.05)
barplot(ego_bp,showCategory = 18,title="The GO_BP enrichment analysis of all DEGs ")+ 
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(ego_bp) str_wrap(ego_bp,width = 25))

KEGG enrichment

library(stringr)
kk<-enrichKEGG(gene      =gene.df$ENTREZID,
               organism = 'hsa',
               pvalueCutoff = 0.05)

> kk[1:30]
               ID                                          Description GeneRatio  BgRatio
hsa04610 hsa04610                  Complement and coagulation cascades   38/1407  79/7431
hsa04933 hsa04933 AGE-RAGE signaling pathway in diabetic complications   40/1407  99/7431
hsa00071 hsa00071                               Fatty acid degradation   23/1407  44/7431
hsa04974 hsa04974                     Protein digestion and absorption   35/1407  90/7431
hsa05146 hsa05146                                           Amoebiasis   36/1407  96/7431
hsa04978 hsa04978                                   Mineral absorption   23/1407  51/7431
hsa04115 hsa04115                                p53 signaling pathway   29/1407  72/7431
hsa05144 hsa05144                                              Malaria   22/1407  49/7431
hsa00982 hsa00982                    Drug metabolism - cytochrome P450   28/1407  72/7431
hsa00980 hsa00980         Metabolism of xenobiotics by cytochrome P450   29/1407  76/7431
hsa03320 hsa03320                               PPAR signaling pathway   28/1407  74/7431
hsa05204 hsa05204                              Chemical carcinogenesis   30/1407  82/7431

mykegg<-barplot(kk,showCategory = 20, title="The KEGG enrichment analysis of all DEGs")+
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(kk) str_wrap(kk,width = 25))
 dotplot(kk,showCategory = 20, title="The KEGG enrichment analysis of all DEGs")+
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(kk) str_wrap(kk,width = 25))

用cowplot拼起来

library(cowplot)
plot_grid(mygobp,mykegg,labels = c("A","B"),align = "h")

Rplot01.jpeg

10 KEGG信号通路的可视化（pathview包）

备注，前面详细部分用gageData重新进行了富集分析，这部分直接调用ClusterProfiler包的富集结果

• 选取enrichKEGG结果pvalue排名第一的hsa4610:Complement and coagulation cascades和排名第7的hsa04115:p53 signaling pathway

library("pathview")
foldchanges = sig.gene$log2FoldChange
names(foldchanges)= gene.df$ENTREZID
head(foldchanges)
pathview(gene.data = foldchanges, pathway.id = "hsa04610", species="hsa")

------------------------------------------------------------------------------

hsa4610:Complement and coagulation cascades

hsa04610.pathview.png

------------------------------------------------------------------------------

hsa04115:p53 signaling pathway

hsa04115.pathview.png

11 把counts结果，clusterprofiler部分的symbol name 和DEGs全部合并到一个表格

备注，这部分的主要问题是没用可以merge用的ID，先看下

> head(mycounts)
                control1 control2 treat1 treat2
ENSG00000000003     4004      781   4229    756
ENSG00000000005        1        0      0      0
ENSG00000000419      776      140   1180    165
ENSG00000000457      624      144   1271    240
ENSG00000000460      260       52    610    105
ENSG00000000938      161       59     57     16

> head(gene.df)
          ENSEMBL   SYMBOL ENTREZID
1 ENSG00000130649   CYP2E1     1571
3 ENSG00000105697     HAMP    57817
4 ENSG00000167244     IGF2     3481
5 ENSG00000162366 PDZK1IP1    10158
6 ENSG00000169715     MT1E     4493
7 ENSG00000074276    CDHR2    54825

> head(sig.gene)
log2 fold change (MLE): condition control vs treat 
Wald test p-value: condition control vs treat 
DataFrame with 6 rows and 6 columns
                  ENSEMBL log2FoldChange     lfcSE      stat        pvalue          padj
                <numeric>      <numeric> <numeric> <numeric>     <numeric>     <numeric>
ENSG00000130649 78059.342       7.190959 0.2110024  34.07998 1.460506e-254 3.017990e-250
ENSG00000268925  5666.055       7.534595 0.2573299  29.27991 1.869559e-188 1.931628e-184
ENSG00000105697  8791.483       6.611356 0.2338215  28.27523 6.970959e-176 4.801597e-172
ENSG00000167244  8988.099       6.767261 0.2582403  26.20528 2.313301e-151 1.195051e-147
ENSG00000162366  4037.960      -7.076854 0.2704122 -26.17062 5.742578e-151 2.373293e-147
ENSG00000169715  3803.597       6.238596 0.2403373  25.95767 1.489813e-148 5.130914e-145

所以执行以下代码：

head(sig.gene)
head(gene.df)
ENSEMBL<-rownames(sig.gene)
sig.gene<-cbind(ENSEMBL, sig.gene)
colnames(sig.gene)[1]<-c("ENSEMBL")

head(mycounts)
ENSEMBL<-rownames(mycounts)
mycounts<-cbind(ENSEMBL, mycounts)
colnames(mycounts)[1]<-c("ENSEMBL")

merge1<-merge(sig.gene,mycounts,by="ENSEMBL")
merge2<-merge(merge1,gene.df, by="ENSEMBL")
head(merge2)
write.csv(DEG_symbole,file="hcc_DEGs_last_results.csv")

结果如下：

> head(merge2)
DataFrame with 6 rows and 13 columns
          ENSEMBL   baseMean log2FoldChange     lfcSE      stat       pvalue         padj  control1  control2    treat1    treat2
      <character>  <numeric>      <numeric> <numeric> <numeric>    <numeric>    <numeric> <integer> <integer> <integer> <integer>
1 ENSG00000000938   70.51956       2.169409 0.4610591  4.705274 2.535248e-06 2.632581e-05       161        59        57        16
2 ENSG00000001460   77.93853      -1.496365 0.4247974 -3.522537 4.274371e-04 2.517116e-03        68        18       268        68
3 ENSG00000001561  382.90350      -1.466940 0.3397462 -4.317751 1.576269e-05 1.347065e-04       416        69      1780       227
4 ENSG00000001626   70.41775       4.713865 0.5812437  8.109965 5.063455e-16 1.944819e-14       226        61        16         1
5 ENSG00000001630  165.86393      -2.283850 0.4094980 -5.577194 2.444286e-08 3.665365e-07        84        28       442       207
6 ENSG00000002549 1683.37557       1.519578 0.2216429  6.855972 7.082930e-12 1.708432e-10      4166      1104      2171       481
       SYMBOL    ENTREZID
  <character> <character>
1         FGR        2268
2       STPG1       90529
3       ENPP4       22875
4        CFTR        1080
5     CYP51A1        1595
6        LAP3       51056

至此，这部分结果可以用来做其他很多下游分析了。