# 学员的转录组实战笔记之胰腺癌差异

Posted on by

文章是：《Identification and functional analysis of novel oncogene DDX60L in pancreatic ductal adenocarcinoma》，细读它可以发现两个不同的数据集， (GSE171485 and GSE171486).

  • SRP313524-PRJNA719796-GSE171485 （polyadenylated-RNA extracted from 6 PDAC specimens and 6 non-tumor adjacent tissues.）
  • SRP313537-PRJNA719795-GSE171486 （knockdown MYEOV, KCNN4, S100A16 and DDX60L in MiaPaCa2 and PANC-1 cells.）

第一个数据集是胰腺癌的癌症和癌旁或者其它对照组织差异，就12个样品，处理起来比较方便，第二个数据集样品数量稍微有一点点多，后面有机会再处理它。文章描述的转录组测序数据的生物信息学处理方法非常陈旧了：

生物信息学处理方法非常陈旧

我们一般来说不会再使用这样的软件和流程，而且作者主要的可视化其数据分析结果也很诡异：

数据分析结果也很诡异

首先呢，我们应该是对这个表达量矩阵进行3张图的质量控制，我在生信技能树的教程：《你确定你的差异基因找对了吗？》提到过，必须要对你的转录水平的全局表达矩阵做好质量控制，最好是看到标准3张图：

  • 左边的热图，说明我们实验的两个分组，正常组织和肿瘤组织的很多基因表达量是有明显差异的
  • 中间的PCA图，说明我们的正常组织和肿瘤组织两个分组非常明显的差异
  • 右边的层次聚类也是如此，说明我们的正常组织和肿瘤组织两个分组非常明显的差异

如果分组在3张图里面体现不出来，实际上后续差异分析是有风险的。这个时候需要根据你自己不合格的3张图，仔细探索哪些样本是离群点，自行查询中间过程可能的问题所在，或者检查是否有其它混杂因素，都是会影响我们的差异分析结果的生物学解释。

但是文章仅仅是拿3个例子说明了，肿瘤样品和肿瘤样品的相关性，正常样品和正常样品的相关性，都是很高的，而且肿瘤和正常相关性很低，这样就是为了说明正常组织和肿瘤组织两个分组的组间差异是大于它们各自内部的组内差异。而且文章里面的差异分析后基因：607 genes were upregulated and 763 genes were downregulated in PDAC specimens using cutoff criteria of |log2(fold- change [FC])| > 0.5 (log2(FC) > 0.5) and p < 0.05

热图显示其实有一些样品是离群点，但是作者并没有做处理。

首先使用conda安装必备软件

conda create -n rna python=3.8
conda activate rna
conda install -y -c bioconda fastqc samtools hisat2 subread 
conda clean -i 
conda install -c "bioconda/label/cf201901" trim-galore

接下来可以使用如下所示代码检验是否安装成功；

fastqc --help 1>/dev/null
trim_galore --help 1>/dev/null
samtools --help 1>/dev/null
hisat2 --help 1>/dev/null
featureCounts --help 2>/dev/null

发现前面安装的trim-galore版本有问题，后面流程报错了，仔细排查后重新确定版本：

conda install -y trim-galore==0.6.7

然后可以使用ascp下载公共数据集

比如这个PRJNA719796的实战，需要自己根据如下链接去EBI里面搜索到，然后自己构建一个 fq.txt 路径文件：

脚本如下：

# conda activate download 
# 自己搭建好 download 这个 conda 的小环境哦。 
# 下面的命令内容构建成为脚本： step1-aspera.sh
cat fq.txt |while read id
do
ascp -QT -l 300m -P33001 \
-i ~/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh \
era-fasp@$id .
done
# nohup bash step1-aspera.sh 1>step1-aspera.log 2>&1 &

这个脚本会根据你在EBI里面搜索到的 fq.txt 路径文件，来批量下载fastq测序数据文件。

查看我们的项目数据（fastq文件）

可以看到， 居然是传说中的单端测序数据，应该是有历史年代感了，如下所示：

ls -lh |cut -d" " -f 5-

6.3G 12月 13 2022 SRR14144009.fastq
5.6G 12月 13 2022 SRR14144010.fastq
5.1G 12月 13 2022 SRR14144011.fastq
6.6G 12月 13 2022 SRR14144012.fastq
6.7G 12月 13 2022 SRR14144013.fastq
6.8G 12月 14 2022 SRR14144014.fastq
3.2G 12月 13 2022 SRR15539279.fastq
2.7G 12月 13 2022 SRR15539280.fastq
4.4G 12月 13 2022 SRR15539281.fastq
3.0G 12月 13 2022 SRR15539282.fastq
2.5G 12月 13 2022 SRR15539283.fastq
5.7G 12月 13 2022 SRR15539284.fastq

首先是质量控制：

代码有的是脚本，有的控制台运行，需要有一点Linux背景才能看懂下面的：

nohup fastqc -t 16 SRR*.fastq &

cat ID | while read id
do
trim_galore -q 20 --length 36 --max_n 3 --stringency 3 --fastqc -o ${cleandata} ${rawdata}/${id}.fastq
done

nohup sh trim_galore.sh >trim_galore.log &

#质控fastp 
cat /home/data/t080304/GEO/SRP313524/3-qc-fastq/ID | while read id 
do
fastp -l 36 -q 20 -W 10 --compression=6 \
-i ${rawdata}/${id}.fastq \
-o ${cleandata}/${id}_clean.fq.gz \
-R ${cleandata}/${id} \
-h ${cleandata}/${id}.fastp.html \
-j ${cleandata}/${id}.fastp.json 
done

nohup sh fastp.sh &

转录组必备数据库文件

主要是下载fa格式的 抽空基因组文件，以及配套的gtf格式的基因组注释文件，如下所示：

# 下载基因组序列
nohup wget -c http://ftp.ensembl.org/pub/release-105/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz &

# 下载转录组序列
nohup wget -c http://ftp.ensembl.org/pub/release-105/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz &
# 下载基因组注释文件
nohup wget -c http://ftp.ensembl.org/pub/release-105/gtf/homo_sapiens/Homo_sapiens.GRCh38.105.chr.gtf.gz &
nohup wget -c http://ftp.ensembl.org/pub/release-105/gff3/homo_sapiens/Homo_sapiens.GRCh38.105.chr.gff3.gz &
# 解压
nohup gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz Homo_sapiens.GRCh38.cdna.all.fa.gz &

# Hisat2构建索引，构建索引时间比较长，建议提交后台运行
nohup hisat2-build Homo_sapiens.GRCh38.dna.primary_assembly.fa GRCh38.dna &

Hisat2构建索引，构建索引时间比较长，建议提交后台运行，而且需要等待起码两个小时以上哈。

使用的是Hisat2比对，而且是单端测序数据的摸索，如下所示：

cat ID | while read id
do
hisat2 -p 16 -x ${index} -U ${inputdir}/${id}_clean.fq.gz -S ${outdir}/${id}.Hisat_aln.sam
samtools sort -@ 15 -o ${outdir}/${id}.Hisat_aln.bam ${outdir}/${id}.Hisat_aln.sam
done

全部的bam文件拿到了之后，进行一次featureCounts对bam文件进行计数

# featureCounts对bam文件进行计数
featureCounts -T 6 -p -t exon -g gene_id -a $gtf -o all.id.txt $inputdir/*_aln.bam

就可以拿到了表达量矩阵，可以完美的进行下游分析啦，读入这个表达量矩阵文件( all.id.txt )

三张图

可以看到并不是完美的区分两个组别，反而是很多混淆的，所以后面的差异分析也很奇葩，如果我们是仅仅是使用默认参数，就会导致根本就没有差异基因（统计学显著的），如下所示：

根本就没有差异基因

这个时候有两个可能性，首先是我们的样品处理阶段，有一些病人取样的时候癌症和癌旁组织可能会混淆，所以需要删除一些样品，6 vs 6 其实是很多可以操作的空间，比如每个组内删除2个样品，这样就是4 vs 4 的组合进行差异分析。

邀请优秀的小伙伴讲解自己如何完成学徒作业

目前学徒作业已经持续了四年多，而且有两三百个不同方向的练习题，我们会邀请小伙伴参与分享自己的完成学徒作业的经历。本次分享的就是这个学员的转录组实战笔记之胰腺癌差异，视频号见：

​

Comments are closed.