CNS图表复现10—表达矩阵是如何得到的

分享是一种态度

前言

CNS图表复现之旅前面我们已经进行了9讲，你可以点击图表复现话题回顾。如果你感兴趣也想加入交流群，自己去：你要的rmarkdown文献图表复现全套代码来了（单细胞）找到我们的拉群小助手哈。

前面的教程里面：CNS图表复现07—原来这篇文章有两个单细胞表达矩阵，我们提到过，是自己读取作者上传到谷歌云里面的2个csv表达矩阵，这个时候有读者就提出来了疑问，作者是如何拿到表达矩阵的呢？

其实呢，这个就涉及到了RNA-seq数据分析的上游流程，需要一些Linux知识啦， 如果你没有服务器，下面的教程就纯粹看一眼哈。

在Linux服务器里面安装conda以及配置aspera下载环境

如果是全新服务器或者全新用户，首先需要安装conda（最适合初学者的软件管理解决方案）：

#一路yes下去
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-4.6.14-Linux-x86_64.sh
source ~/.bashrc

然后使用conda安装一些软件或者软件环境，比如下载测序数据文件的aspera软件环境：

conda create -n download -y 
conda activate download 
conda install -y -c hcc aspera-cli 
which ascp 
## 一定要搞清楚你的软件被conda安装在哪
ls -lh ~/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh

根据文章找到数据下载链接

参考：使用ebi数据库直接下载fastq测序数据  , 需要自行配置好，然后去EBI里面搜索到的 fq.txt 路径文件：

• 项目地址是 ：https://www.ebi.ac.uk/ena/browser/view/PRJNA591860

可以看到是超过2万个文件哦：

拿到全部的下载链接，存储为文本文件（ fq.txt），如下：

fasp.sra.ebi.ac.uk:/vol1/fastq/SRR107/015/SRR10777215/SRR10777215_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR107/015/SRR10777215/SRR10777215_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR107/016/SRR10777216/SRR10777216_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR107/016/SRR10777216/SRR10777216_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR107/017/SRR10777217/SRR10777217_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR107/017/SRR10777217/SRR10777217_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR107/018/SRR10777218/SRR10777218_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR107/018/SRR10777218/SRR10777218_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR107/019/SRR10777219/SRR10777219_1.fastq.gz

然后使用简单的脚本，读取那个文本文件（ fq.txt）进行批量下载，脚本如下：

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 &

可以看到下载的文件如下：

 194M Oct  4 11:07 raw/SRR10777215_1.fastq.gz
 200M Oct  4 11:09 raw/SRR10777215_2.fastq.gz
  92M Oct  4 11:09 raw/SRR10777216_1.fastq.gz
  94M Oct  4 11:10 raw/SRR10777216_2.fastq.gz
  40M Oct  4 11:11 raw/SRR10777217_1.fastq.gz
  41M Oct  4 11:12 raw/SRR10777217_2.fastq.gz
  51M Oct  4 11:13 raw/SRR10777218_1.fastq.gz
  
  ## 中间省略一万多个文件：
  
  52M Oct 13 09:55 raw/SRR10783212_2.fastq.gz
  51M Oct 13 09:55 raw/SRR10783212_1.fastq.gz
 120M Oct 13 09:54 raw/SRR10783211_2.fastq.gz
 116M Oct 13 09:53 raw/SRR10783211_1.fastq.gz
 175M Oct 13 09:52 raw/SRR10783210_2.fastq.gz
 161M Oct 13 09:51 raw/SRR10783210_1.fastq.gz
  15M Oct 13 09:50 raw/SRR10783209_2.fastq.gz
  16M Oct 13 09:49 raw/SRR10783209_1.fastq.gz

横跨了一个国庆节假期，才下载了不到一半的文件。因为仅仅是演示这个项目的流程，所以我就直接终止这个程序啦。如果上面的代码你完全看不懂呢，需要自己去看B站视频，详见：我在b站的74小时生信工程师教学视频合辑：

都是需要R语言和linux基础的哦！

然后走最简单的hisat2+featureCounts流程拿到表达矩阵

我这里简单的演示几个单细胞转录组样本即可，都是双端测序，如下所示：

$ls -lh  raw/SRR1077721* |cut -d" " -f 5-
194M Oct  4 11:07 raw/SRR10777215_1.fastq.gz
200M Oct  4 11:09 raw/SRR10777215_2.fastq.gz
 92M Oct  4 11:09 raw/SRR10777216_1.fastq.gz
 94M Oct  4 11:10 raw/SRR10777216_2.fastq.gz
 40M Oct  4 11:11 raw/SRR10777217_1.fastq.gz
 41M Oct  4 11:12 raw/SRR10777217_2.fastq.gz
 51M Oct  4 11:13 raw/SRR10777218_1.fastq.gz
 54M Oct  4 11:13 raw/SRR10777218_2.fastq.gz
 89M Oct  4 11:16 raw/SRR10777219_1.fastq.gz
 91M Oct  4 11:17 raw/SRR10777219_2.fastq.gz

一般来说raw的fq文件，需要走一下trim_galore软件进行质控：

conda activate qc
mkdir -p clean
trim_galore --help

for id in  {15..19} 
do 
    nohup  trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired \
    -o  clean   raw/SRR107772$id*.fastq.gz   & 
done 

得到的干净的fq文件如下：

189M Oct 13 11:01 clean/SRR10777215_1_val_1.fq.gz
194M Oct 13 11:01 clean/SRR10777215_2_val_2.fq.gz
 90M Oct 13 10:59 clean/SRR10777216_1_val_1.fq.gz
 92M Oct 13 10:59 clean/SRR10777216_2_val_2.fq.gz
 39M Oct 13 10:57 clean/SRR10777217_1_val_1.fq.gz
 40M Oct 13 10:57 clean/SRR10777217_2_val_2.fq.gz
 51M Oct 13 10:58 clean/SRR10777218_1_val_1.fq.gz
 53M Oct 13 10:58 clean/SRR10777218_2_val_2.fq.gz
 82M Oct 13 10:59 clean/SRR10777219_1_val_1.fq.gz
 84M Oct 13 10:59 clean/SRR10777219_2_val_2.fq.gz

可以看到过滤前后的fq文件大小是有变化的。

比对呢，直接选择hisat2软件，需要自己搞定参考基因组索引文件哦，代码如下；

conda activate rna
mkdir -p align
index=$HOME/reference/index/hisat/hg38/genome
for id in  {15..19} 
do 
fq1=clean/SRR107772${id}_1_val_1.fq.gz
fq2=clean/SRR107772${id}_2_val_2.fq.gz
hisat2 -p 4 -x $index -1  $fq1 -2  $fq2 | samtools sort -@ 4  -o align/$id.bam -
done  

得到的bam文件如下：

777M Oct 13 11:24 15.bam
312M Oct 13 11:25 16.bam
137M Oct 13 11:26 17.bam
132M Oct 13 11:26 18.bam
169M Oct 13 11:27 19.bam

可以看到，上面的脚本输出bam文件的时候，我忘记把 SRR107772 这个前缀加上去了。

因为是人类癌症病人数据，featureCounts计数代码如下：

mkdir counts 
cd counts

wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_35/gencode.v35.chr_patch_hapl_scaff.annotation.gtf.gz
gunzip gencode.v35.chr_patch_hapl_scaff.annotation.gtf.gz

gtf=gencode.v35.chr_patch_hapl_scaff.annotation.gtf
 featureCounts -T 4 -p -t exon -g gene_name  -a $gtf -o all.counts.id.txt  \
  *.bam 1>counts.id.log 2>&1 &

这样就拿到了表达矩阵，当然了，如果是全部的几万个文件的运行，耗时就很可观了，耗费的计算资源也不小哦。

当然了，这个代码有点超纲了，对绝大部分初学者来说。需要把Linux的6个阶段跨越过去 ，一般来说，每个阶段都需要至少一天以上的学习：

• 第1阶段：把linux系统玩得跟Windows或者MacOS那样的桌面操作系统一样顺畅，主要目的就是去可视化，熟悉黑白命令行界面，可以仅仅以键盘交互模式完成常规文件夹及文件管理工作。
• 第2阶段：做到文本文件的表格化处理，类似于以键盘交互模式完成Excel表格的排序、计数、筛选、去冗余，查找，切割，替换，合并，补齐，熟练掌握awk,sed,grep这文本处理的三驾马车。
• 第3阶段：元字符，通配符及shell中的各种扩展，从此linux操作不再神秘！
• 第4阶段：高级目录管理：软硬链接，绝对路径和相对路径，环境变量。
• 第5阶段：任务提交及批处理，脚本编写解放你的双手。
• 第6阶段：软件安装及conda管理，让linux系统实用性放飞自我。

详见：《生信分析人员如何系统入门Linux(2019更新版)》

往期回顾

---

---

---

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

看完记得顺手点个“在看”哦！

生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注