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

2020-11-02 10:17:49 浏览数 (1)

分享是一种态度

前言

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

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

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

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

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

代码语言：javascript复制
#一路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软件环境：

代码语言：javascript复制
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），如下：

代码语言：javascript复制
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）进行批量下载，脚本如下：

代码语言：javascript复制
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 &

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

代码语言：javascript复制
 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小时生信工程师教学视频合辑：

  • 免费视频课程《RNA-seq数据分析》
  • 免费视频课程《WES数据分析》
  • 免费视频课程《ChIP-seq数据分析》
  • 免费视频课程《ATAC-seq数据分析》
  • 免费视频课程《TCGA数据库分析实战》
  • 免费视频课程《甲基化芯片数据分析》
  • 免费视频课程《影像组学教学》
  • 免费视频课程《LncRNA-seq数据》
  • 免费视频课程《GEO数据挖掘》

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

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

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

代码语言：javascript复制
$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软件进行质控：

代码语言：javascript复制
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文件如下：

代码语言：javascript复制
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软件，需要自己搞定参考基因组索引文件哦，代码如下；

代码语言：javascript复制
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文件如下：

代码语言：javascript复制
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计数代码如下：

代码语言：javascript复制
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更新版)》

往期回顾

IRIS3 ：细胞类型特异的调节子分析框架

剑桥大学5天功能基因组培训课程全套资料（含单细胞转录组）

二十年前做科研你只需要检测一些基因在一些癌症细胞系表达量情况即可

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

  • 单细胞初级8讲和高级分析8讲
  • 生信爆款入门-第9期（线上直播4周，马拉松式陪伴，带你入门）
  • 数据挖掘第7期（线上直播3周，马拉松式陪伴，带你入门）
数据分析 linux 数据挖掘

0 人点赞