ixxmu
commented
1 year ago 学徒抽丝剥茧想搞清楚这个转录组数据问题出在哪里 by 生信技能树
前些天布置了学徒作业- CNP0002454的分析，详见：
这个CNP0002454数据集来源的文章：Identification of a 3-Gene Model as Prognostic Biomarker in Patients With Gastric Cancer ，大家可以自行阅读
DOI：https://doi.org/10.3389/fonc.2022.930586
0.下载数据
数据在CNGBdb，https://db.cngb.org/search/project/CNP0002454/
查了一下批量下载的教程，发现一般是用Aspera下载
Aspera批量下载：http://www.bio-info-trainee.com/8587.html
但是需要自己手动制作new_fq.txt文件，也就是把下载链接一条条复制粘贴进去。
看了一下本文的样本有34个，每个样本双端测序，也就是要复制粘贴68次，还是有点繁琐了。
本来尝试了另一种傻瓜式的下载方法：
下载Xftp或FileZilla 进入ftp://ftp.cngb.org，匿名登录 通过文件夹路径找到项目文件夹：/pub/CNSA/data4/CNP0002454 把整个CNP0002454文件夹拖进服务器目标路径即可
下载到一半连接断开了……
缺点：下载速度大概3M/s，有点慢，网容易断。
最后，我是用以下方法下载的：
只需要知道项目编号，找到ftp路径：ftp://ftp.cngb.org/pub/CNSA/data4/CNP0002454
然后通过wget命令递归下载，速度20M/s，很快：
wget -r -nH -nd -P ./ ftp://ftp.cngb.org/pub/CNSA/data4/CNP0002454 --ftp-user=anonymous --ftp-password=anonymous@example.com
# -r :递归下载
# -nH：不创建主机目录
# -nd：不创建目录
# -P：将文件保存到目录
这样可以把所有文件都下载在一个文件夹里。
注意：一定要加-nd参数！否则会得到一个超级无敌长的文件夹套娃！
1.质控过滤
FastQC质控，Trim_galore过滤。
首先用FastQC检测原始数据质量：
fastqc -t 12 -o ./ *.gz
# -t：线程数
multiqc ./ -n "All_Raw_fq"
再对原始数据进行过滤及质量再检测。
改了一下代码，在曾老师代码的基础上增加了脚本的泛用性：
# Trim_FastQC.sh
#!/bin/bash
######################
# usage: bash Trim_FastQC.sh <QCconfig> <Output_dir> <Raw_file_dir>
# QCconfig：遍历列表
# Output_dir：输出文件路径
# Raw_file_dir：原始数据文件路径
######################
echo =================== QC Start ===================
echo QC method: Trim Galore + FastQC
# 创建<Output_dir>输出文件夹
if [ -d "$2" ] ; then
echo $2 have been existed.
else
mkdir $2
echo Make $2
fi
cd $2
source /home/miniconda3/bin/activate YAP-ChIP
# 根据QCconfig遍历
cat $1 | while read id;
do
arr=($id)
path1=${arr[0]}
path2=${arr[1]}
fq1=$(basename $path1)
fq2=$(basename $path2)
echo ------------------ ${fq1%%.*} ------------------
echo fq1: $fq1
echo fq2: $fq2
trim_galore --fastqc --gzip --basename ${fq1%%.*} -O $2 --paired $3/$path1 $3/$path2
# 根据文件名可以修改--basename参数
if [ $? -ne 0 ] ; then
echo ------------------ Failed ------------------
else
echo ------------------ Success ------------------
fi
done
# multiqc 查看多样本QC结果
multiqc ./ -n "All_clean_fq"
if [ $? -ne 0 ] ; then
echo ----- MultiQC Failed -----
else
echo ----- MultiQC Success -----
echo ==================== QC End ====================
使用脚本前需要做一个QCconfig文件：
ls *.gz | grep _1.fq.gz > ../fq1
ls *.gz | grep _2.fq.gz > ../fq2
paste fq1 fq2 > QCconfig
# QCconfig
DP8400018285BL_L01_100_1.fq.gz DP8400018285BL_L01_100_2.fq.gz
DP8400018285BL_L01_101_1.fq.gz DP8400018285BL_L01_101_2.fq.gz
DP8400018285BL_L01_102_1.fq.gz DP8400018285BL_L01_102_2.fq.gz
DP8400018285BL_L01_103_1.fq.gz DP8400018285BL_L01_103_2.fq.gz
...
然后运行：
nohup bash Trim_FastQC.sh /data2/RNA/QCconfig /data2/RNA/1.Trim+QC /data2/RNA/RawData &
QC结果
原始数据：
原始数据质量就已经合格了，而且也没有检测到接头污染。
过滤数据：
但是发现GC含量和重复序列水平不太正常。
原始数据：
过滤数据：
只有DP8400018285BL_L01_128勉强通过GC含量和重复序列水平检测。考虑到样品是胃组织，可能存在细菌污染情况。
2.Hisat2比对
比对使用脚本：
# hisat2.sh
#!/bin/bash
# usage: hisat2.sh <Alignconfig> <Rawdata_dirpath> <output_dirpath>
# mm10=/data2/Ref/hisat2_index/mm10/grcm38.p6/GRCm38.p6.genome
hg19=/data2/Ref/hisat2_index/hg19/hg19/genome
echo ========== RNA-seq Alignment ==========
echo Method:Hisat2
if [ -d "$3" ];
then
echo $3 have been existed.
else
mkdir $3
echo Make $3.
fi
cd $3
source /home/miniconda3/bin/activate RNA-seq
cat $1 | while read id;
do
arr=($id)
fq1=${arr[0]}
fq2=${arr[1]}
out=${fq1%%_1_val*}
echo fq1:$fq1
echo fq2:$fq2
hisat2 -x $hg19 -p 20 -1 $2/$fq1 -2 $2/$fq2 | samtools view -Sbh - > $3/${out}.bam
if [ $? -ne 0 ];
then
echo ------------------ Failed ------------------
else
echo ------------------ Success ------------------
fi
done
echo ========== RNA-seq Alignment End ==========
运行之前先制作Alignconfig：
# 在Trim之后的文件夹中
ls *.gz | grep _1.fq.gz > ../fq1
ls *.gz | grep _2.fq.gz > ../fq2
paste fq1 fq2 > Alignconfig
# Alignconfig
DP8400018285BL_L01_100_1_val_1.fq.gz DP8400018285BL_L01_100_1_val_2.fq.gz
DP8400018285BL_L01_101_1_val_1.fq.gz DP8400018285BL_L01_101_1_val_2.fq.gz
DP8400018285BL_L01_102_1_val_1.fq.gz DP8400018285BL_L01_102_1_val_2.fq.gz
DP8400018285BL_L01_103_1_val_1.fq.gz DP8400018285BL_L01_103_1_val_2.fq.gz
...
然后运行hisat2.sh：
nohup bash hisat2.sh /data2/Alignconfig /data2/raw /data2/2.Hisat2 & 1>hisat2.log
检查比对率：
# 图省事可以一行命令
cat hisat2.log | grep 'fq1\|overall alignment rate' | sed 's/fq1:\|_1_.*\|over.*//g' > AlignRate.txt
# 分两步方便查看
cat hisat2.log | grep fq1 | sed 's/fq1:\|_1_.*//g' > ID
cat hisat2.log | grep overall| sed 's/over.*//g' > Alignrate
paste ID Alignrate > AlignRate.txt
# 按比对率大小排序
sort -n -k 2 -r AlignRate.txt
可以发现大多数样本的比对率在70-80%，个别低至46%，因为前面QC看到了大量非人类物种的序列干扰，所以我们使用人类参考基因组比对那些reads肯定是会失败。（正常情况下我们的转录组数据人或者小鼠比对都是99%以上）
而GC含量勉强合格的DP8400018285BL_L01_128比对率86.42%是最高的。
# AlignRate.txt
DP8400018285BL_L01_128 86.42%
DP8400018285BL_L01_125 84.55%
DP8400018285BL_L01_117 83.53%
DP8400018285BL_L01_98 82.01%
DP8400018285BL_L01_127 81.99%
DP8400018285BL_L01_26 81.75%
DP8400018285BL_L01_91 81.29%
DP8400018285BL_L01_100 81.15%
DP8400018285BL_L01_103 80.87%
DP8400018285BL_L01_104 80.79%
DP8400018285BL_L01_99 80.77%
DP8400018285BL_L01_123 80.76%
DP8400018285BL_L01_102 80.36%
DP8400018285BL_L01_85 80.03%
DP8400018285BL_L01_25 79.94%
DP8400018285BL_L01_101 79.91%
DP8400018285BL_L01_30 79.74%
DP8400018285BL_L01_29 79.17%
DP8400018285BL_L01_96 79.06%
DP8400018285BL_L01_124 78.86%
DP8400018285BL_L01_92 78.16%
DP8400018285BL_L01_97 75.84%
DP8400018285BL_L01_94 75.72%
DP8400018285BL_L01_86 75.68%
DP8400018285BL_L01_126 75.08%
DP8400018285BL_L01_28 75.01%
DP8400018285BL_L01_87 73.48%
DP8400018285BL_L01_122 73.13%
DP8400018285BL_L01_121 72.98%
DP8400018285BL_L01_88 72.05%
DP8400018285BL_L01_93 68.29%
DP8400018285BL_L01_90 63.78%
DP8400018285BL_L01_89 57.57%
DP8400018285BL_L01_95 46.37%
要进行污染检测可以使用Blast或者FastQ Screen。
Blast需要下载庞大的数据库文件，FastQ Screen需要下载bowtie、bowtie2或者bwa的物种索引文件。
需要的数据库文件太大了这里就不做了，放一些攻略：
Blast:
https://www.jianshu.com/p/2487cfe27e0c https://blog.csdn.net/qq_42962326/article/details/105081327
FastQ Screen:
https://www.jianshu.com/p/0990f478ef63
3.featureCounts定量
计算比对过的文件中位于exon部分的reads数。这个部分应该可以过滤掉大部分污染，只计算能够映射到转录本的reads。
featureCounts -T 12 -p -t exon -g gene_id \
-a ../Annotation/GRCh37.p13/gencode.v19.annotation.gtf \
--extraAttributes gene_name -o RNA.txt \
*.bam \
2>featureCounts.log
Tip：加--extraAttributes参数把gene_name(也就是Symbol)一起输出，在后面在R中就不用转换了。（其实 -g 参数也可以指定 使用gene_name(也就是Symbol) ）
运行结束后对映射结果进行质控：
multiqc ./ -n "hisat2QC"
得到映射率：
DP8400018285BL_L01_100 40.0%
DP8400018285BL_L01_101 34.7%
DP8400018285BL_L01_102 34.1%
DP8400018285BL_L01_103 37.3%
DP8400018285BL_L01_104 36.6%
DP8400018285BL_L01_117 39.3%
DP8400018285BL_L01_121 32.9%
DP8400018285BL_L01_122 34.2%
DP8400018285BL_L01_123 39.1%
DP8400018285BL_L01_124 35.3%
DP8400018285BL_L01_125 26.0%
DP8400018285BL_L01_126 35.1%
DP8400018285BL_L01_127 40.9%
DP8400018285BL_L01_128 36.8%
DP8400018285BL_L01_25 43.4%
DP8400018285BL_L01_26 39.0%
DP8400018285BL_L01_28 42.2%
DP8400018285BL_L01_29 44.3%
DP8400018285BL_L01_30 37.6%
DP8400018285BL_L01_85 20.7%
DP8400018285BL_L01_86 20.8%
DP8400018285BL_L01_87 23.2%
DP8400018285BL_L01_88 19.0%
DP8400018285BL_L01_89 13.1%
DP8400018285BL_L01_90 14.9%
DP8400018285BL_L01_91 37.3%
DP8400018285BL_L01_92 27.8%
DP8400018285BL_L01_93 24.1%
DP8400018285BL_L01_94 26.2%
DP8400018285BL_L01_95 18.7%
DP8400018285BL_L01_96 19.4%
DP8400018285BL_L01_97 22.4%
DP8400018285BL_L01_98 26.2%
DP8400018285BL_L01_99 41.0%
可以看到这个数据集不仅仅是前面的参考基因组比对率非常低，而且呢，它映射到参考基因组的基因注释区域比例也很多，正常情况下我们的转录组可以达到60%~80%的映射率。
另外，由于做差异分析之前要生成一个group_list标记样本的处理对照组关系，在featureCounts运行结束之后把样本id改成样本类型。
先做一个id_sample.txt对应表：
# id_sample.txt
DP8400018285BL_L01_100.bam case1
DP8400018285BL_L01_101.bam control1
DP8400018285BL_L01_102.bam case2
DP8400018285BL_L01_103.bam control2
DP8400018285BL_L01_104.bam case3
DP8400018285BL_L01_117.bam control3
DP8400018285BL_L01_121.bam control4
DP8400018285BL_L01_122.bam case4
DP8400018285BL_L01_123.bam control5
DP8400018285BL_L01_124.bam case5
DP8400018285BL_L01_125.bam control6
DP8400018285BL_L01_126.bam case6
...
再写一个循环即可：
cat id_sample.txt | while read id;
do
arr=($id)
id=${arr[0]}
sample=${arr[1]}
# -i表示改动原文件，加后缀表示修改原文件的同时生成一个拥有该后缀的备份
sed "s/$id/$sample/g" RNA.txt -i.bak
done
4.DESeq2差异分析
以下代码都在R中跑。
首先走DESeq2差异分析流程，得到差异分析结果和DESeq2的标准化表达矩阵：
rm(list = ls())
options(stringsAsFactors = F)
library(openxlsx)
library(patchwork)
library(stringr)
# 去掉ensembl_id中.及之后的字符
a=read.table("../RNA.txt", header = T)
a$Geneid <- str_split(a$Geneid,'\\.',simplify = T)[,1]
exprSetB <- a[,-c(2:6)]
if(F){
library(org.Hs.eg.db)
colnames(exprSetB)[1] <- "ensembl_id"
table(table(exprSetB$ensembl_id)>1)
table(exprSetB$ensembl_id)[table(exprSetB$ensembl_id)>1] #查找重复出现的基因名
exprSetB=exprSetB[order(exprSetB$ensembl_id),] #按ensemblID重新排序
exprSetB=exprSetB[!duplicated(exprSetB$ensembl_id),] #ensembl_id去重
rownames(exprSetB) <- exprSetB$ensembl_id
}
exprSet <- exprSetB
rownames(exprSet) <- exprSetB$ensembl_id
exprSet <- exprSet[,-c(1,2)]
exprSet <- exprSet[,order(colnames(exprSet))] # 列排序
# 差异分析
suppressMessages(library(DESeq2))
group_list <- factor(c(
rep('Case',17),
rep('Control',17)
))
colData <- data.frame(row.names=colnames(exprSet),
group_list=group_list
)
colData$group_list <- relevel(colData$group_list, ref = "Control")
dds <- DESeqDataSetFromMatrix(countData = exprSet,
colData = colData,
design = ~ group_list)
dds <- DESeq(dds)
normalized_counts<-counts(dds,normalized=TRUE) #返回标准化的值
normalized_counts <- as.data.frame(normalized_counts)
write.xlsx(normalized_counts, file = '../normalized_counts.xlsx', rowNames = TRUE)
# 返回差异分析结果
res <- results(dds)
resOrdered <- res[order(res$padj),]
resOrdered <- na.omit(as.data.frame(resOrdered))
write.xlsx(resOrdered,"../DEG_result.xlsx", rowNames = TRUE)
可视化部分
1.PCA图
library(ggplot2)
library(factoextra)
# 准备exprSet表达矩阵
exprSet1 <- t(normalized_counts)
exprSet1 <- exprSet1[,colSums(exprSet1)!=0]
expr_pca <- prcomp(exprSet1, center = T, scale. = T)
expr_pcasum <- summary(expr_pca)
pic_expr_pca <- fviz_pca_ind(expr_pca,
col.ind=group_list,
repel = TRUE, # 避免标签重叠
addEllipses = T,
legend.title="Groups",
ellipse.type="confidence",
#ellipse.level=0.9,
palette = c("#8665C8", "#FF8635"),
title = 'PCA for case&Control'
)+
theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"))#加个边框
ggsave(pic_expr_pca,filename = "../PCA.pdf",width = 6,height = 4)
看PCA的结果感觉差异不是很明显……
PS：其实这个时候，我们应该是做三张图~
我在生信技能树的教程：《你确定你的差异基因找对了吗？》提到过，必须要对你的转录水平的全局表达矩阵做好质量控制，最好是看到标准3张图：
左边的热图，说明我们实验的两个分组，normal和npc的很多基因表达量是有明显差异的 中间的PCA图，说明我们的normal和npc两个分组非常明显的差异 右边的层次聚类也是如此，说明我们的normal和npc两个分组非常明显的差异
如果分组在3张图里面体现不出来，实际上后续差异分析是有风险的。这个时候需要根据你自己不合格的3张图，仔细探索哪些样本是离群点，自行查询中间过程可能的问题所在，或者检查是否有其它混杂因素，都是会影响我们的差异分析结果的生物学解释。
而这个CNP0002454数据集很明显如果是使用case和control这样的分组，三张图是失败的，其实这个时候可以加入我们前面的比对率和映射率，就可以探索为什么会失败。
2.火山图
# plot
nrDEG=resOrdered
logFC_cutoff=1
nrDEG$change = as.factor(ifelse(nrDEG$padj < 0.05 & abs(nrDEG$log2FoldChange) > logFC_cutoff,
ifelse(nrDEG$log2FoldChange > logFC_cutoff ,'UP','DOWN'),'NOT'))
title_vol <- paste0('\nCutoff for log2FC is ',logFC_cutoff,
'\nThe number of up gene is ',nrow(nrDEG[nrDEG$change =='UP',]) ,
'\nThe number of down gene is ',nrow(nrDEG[nrDEG$change =='DOWN',]))
# enhancedvolcano作图
library(EnhancedVolcano)
# 设置作图参数
keyvals = list()
keycol <- function(x){
ifelse(x !='NOT', ifelse(x =='UP', 'red2', 'royalblue'),'grey')
}
keyvals <- lapply(nrDEG$change, keycol)
names(keyvals)[keyvals == 'red2'] <- 'UP'
names(keyvals)[keyvals =='grey' ] <- 'NOT'
names(keyvals)[keyvals == 'royalblue'] <- 'DOWN'
# 火山图
volc_plot <- EnhancedVolcano(nrDEG,
max.overlaps = 10,
lab = rownames(nrDEG),
x = 'log2FoldChange',
y = 'padj',
selectLab = '', # 标记基因
##设置题目
title = "Volcano",
titleLabSize = 20,
subtitle = title_vol,
subtitleLabSize = 16,
legendPosition = 'right',
caption = '',
#captionLabSize = 14,
##设置p值cutoff和foldchange的cutoff
pCutoff = 0.05, FCcutoff = logFC_cutoff,
##设置点和标签的大小
pointSize = 2.0, labSize = 6.0,
#添加Connectors
# drawConnectors = TRUE, widthConnectors = 0.5,
# 四象限颜色
# col = c('grey', 'grey', 'grey', 'red2'),
# 自定义颜色
colCustom = keyvals)
ggsave(volc_plot,filename = "../volcano.pdf",width = 8,height = 6)
分析得到差异基因共557个，其中，上调193个，下调364个。
3.表达量热图
## heatmap
library(pheatmap)
#选择上下调基因
choose_gene <- subset(nrDEG,nrDEG$change != 'NOT')
choose_matrix=normalized_counts[rownames(choose_gene)[order(choose_gene$change, decreasing = TRUE)],]
# 添加注释条
Group <- data.frame(Group=group_list, row.names = colnames(normalized_counts))
Up_Down <- data.frame(Up_Down = choose_gene$change, row.names = rownames(choose_gene))
heat_plot <- pheatmap(choose_matrix,
scale = "row",
cluster_rows = F,
cluster_cols = F,
show_colnames =F,
show_rownames = F,
annotation_names_row = F,
annotation_names_col = F,
annotation_col=Group,
annotation_colors = list(Group = c(Control = "#FFBE7A", Case = "#BEB8DC"),
Up_Down = c(UP = "#F84B0C", DOWN = "#0C79B8")),
annotation_row = Up_Down,
border_color = NULL,
main = 'Pheatmap 557 DEGs',
#angle_col = 315,
color = colorRampPalette(c("#0C79B8", "#FFFCCE", "#F84B0C"))(100))
ggsave(filename = "../heatmap.pdf", plot = heat_plot, width = 5, height = 5)
复现图（上）及原图（下）
虽然差异基因个数有差别，但表达量趋势是相似的。其实这个时候基本上可以判断出来 华大基因单细胞团队的这个差异分析后的热图真奇怪，是因为这个转录组测序数据质量差的问题，比对率和映射率都不好，所以表达量矩阵就有问题，那么后续强行找差异后的可视化也是不对劲。
4.韦恩图
文章分析得到321个差异基因，可以在补充表3中获得所有差异基因列表，基于此可以作韦恩图看两个差异基因集的关系：
library(openxlsx)
library(ggplot2)
library(ggvenn)
# 导入DEGs表
a0 <- read.xlsx(xlsxFile = '../Data Sheet 2.xlsx', sheet = 4,
colNames = TRUE, rowNames = TRUE)
colnames(a0) <- a0[1,]
a0 <- a0[-1,]
# 提取DEGs
a0_DEGs <- rownames(a0)
a_DEGs <- exprSetB[rownames(choose_gene),"gene_name"]
# 作Venn图
Venn <- ggvenn(list(DEGs_321 = a0_DEGs, DEGs_557 = a_DEGs),
fill_color = c("#957BCA", "#FA5D45"), fill_alpha = 0.7,
stroke_color = 'black')
ggsave('../Venn.pdf', plot = Venn, width = 5, height = 5)
原文（左，321个）及复现（右，557个）
交集部分相对于文章找到的差异基因集来说还是很大的。当然了，两次差异分析（我们的复现和文章自己的差异）的对比其实不仅仅是这样的韦恩图看交集。
初学者最喜欢怀疑自己的分析是否正确，比如差异分析的时候就容易陷入上下调基因数量的对比问题，文章可能是上下调一千附近，但是学员自己复现的时候就数量对不上。
其实这个问题并不在于上下调基因数量，应该是看质量，这样的对比才有意义。详见：两次差异分析结果的比较不要局限于韦恩图
思考
不知道污染对后续的分析有没有影响（虽然在featureCounts的时候只计算映射到exon部分的reads），可以污染检测之后用FastQ Screen或者其他软件把污染物种reads删除，再走一次转录组流程对比结果； 比文章得到的DEGs多两百多个，可以继续做GO和GSEA-KEGG分析，与文章结果对比。
https://mp.weixin.qq.com/s/ly2p0gS2wvU7V8s0zeSMHA