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

Posted on by

前面的教程里面：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 路径文件：

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

image-20201013095344766

拿到全部的下载链接，存储为文本文件（ 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

## 中间省略4万多个文件：

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更新版)》

​

Comments are closed.