看优秀本科生如何一周内学会Linux进而搞定RNA-seq上游分析

生信技能树

发布于 2020-04-14 12:46:57

6.5K1

发布于 2020-04-14 12:46:57

距离公布要带500个优秀本科生入门生物信息学的活动不到一个月，虽然真正入选不到一百，但是培养成绩喜人，出勤率接近百分之百，大部分人在短短两个星期就完成了R基础知识学习，Linux认知，甚至看完了转录组实战水平，进而完成了一个自己的课题！如果你也感兴趣这个活动，那么，直达文末找到活动链接，免费申请加入吧！

优秀本科生

谢天 2020-04-07

我是武汉大学基础医学专业第一届的学生，2016年9月刚进大学的时候就选了导师进入实验室接受科研训练。虽然我们实验室不是专门做生物信息学的，但第一次和导师正式交流的时候，她就建议我要学点生信。（巧合的是2016年9月也是生信菜鸟团转型生信技能树的时间点，如果所有的导师都如此明智就好了）

刚开始主要是跟着实验室师兄师姐学了些GEO数据库的芯片数据分析，主要是利用Excel、统计软件以及在线工具分析。慢慢地，我发现只会这些是远远不够的，于是就关注了生信技能树、生信菜鸟团等公众号，希望能学习一些更专业的技能。（比如R语言编程）

我们实验室主要是需要基于GEO或TCGA数据库的转录组数据做一些下游分析，利用现成的工具或R语言分析一下，算不上太难。这两年半零零碎碎地学了点R语言，真正要实操的时候也还是会有诸多问题。

从今年开始，课就比较少了，本来是计划好好训练一下湿实验的技能，但是因为疫情暂时不能实现。所以疫情期间的课外学习计划调整为生信与编程（R、Python），统计与临床研究（包括meta分析），生物学知识补习（组学、前沿 三大学习板块。

正好生信技能树Jimmy老师做了一个本科生毕业设计的辅导活动，虽然我是五年制，明年才毕业，但正好学一下专业的生信分析也在计划之内，所以就报名了。感谢Jimmy老师提供这个宝贵的学习机会?

本来自己暂时没有计划学习Linux（因为确实畏难），但是这段时间跟着jimmy老师教学团队的学习让我明白了为什么Linux是生信分析的基础。以前我们做转录组分析主要是基于表达量的下游分析，对于RNA-seq上游分析几乎是零基础，在听了老师的课之后，自己摸爬滚打了一周，终于勉强走了一遍。下面分享一下我的坎坷历程：

RNA-seq数据GSE140739上游分析记录

周五下午从学校超算中心申请到了一个账号，终于可以在服务器上实操一下了。这次实战用的数据为刘胡丹老师二月份Cancer Cell文章的RNA-seq数据GSE140739。

一、原始测序数据

00 软件安装

软件安装没有什么太特别的地方，按照网上的教程一步一步没有太大问题。?conda管理生信软件一文就够

这里摘录的方法来自《原创10000+生信教程大神给你的RNA实战视频演练》。

1. 安装conda

推荐使用偷懒方法，比如安装miniconda软件，下载地址：https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/ 这样就可以使用它安装绝大部分其它软件。

但是在中国大陆的小伙伴，需要更改镜像源配置。

conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes

2. 安装软件

为了避免污染Linux工作环境，推荐在conda中创建各个流程的安装环境，比如：

conda create -n rna python=3 #创建名为rna的软件安装环境
conda info --envs #查看当前conda环境
source activate rna #激活conda的rna环境

转录组分析需要用到的软件列表

质控 fastqc , multiqc, trimmomatic, cutadapt, trim-galore 比对 star, hisat2, bowtie2, tophat, bwa, subread 计数 htseq, bedtools, deeptools, salmon

安装软件

conda install -y sra-tools
conda install -y trimmomatic
conda install -y cutadapt multiqc 
conda install -y trim-galore
conda install -y star hisat2 bowtie2
conda install -y subread tophat htseq bedtools deeptools
conda install -y salmon
source deactivate #注销当前的rna环境

01 数据下载

1. 下载SRA数据

数据下载第一步就遇到了很多障碍，利用aspera下载的时候遇到了这个报错，在网上也没找到很好的解决方法，从周五晚上一直试到周六晚饭前

ascp: Failed to open TCP connection for SSH, exiting.  
Session Stop  (Error: Failed to open TCP connection for SSH)

在准备放弃的时候，看到了生信技能树的这篇文章《原创10000+生信教程大神给你的RNA实战视频演练》，用prefetch成功下载SRA数据。

从PRJNA590732下载到SRR_Acc_List.txt的文档，文档内容是SRR号码（一个号码占据一行），存放到/dat01/xietian/sra路径下。

cat SRR_Acc_List.txt | while read id; do (prefetch  ${id} );done

下载得到的文件：

├── SRR10502962.sra
├── SRR10502963.sra
├── SRR10502964.sra
├── SRR10502965.sra
├── SRR10502966.sra
└── SRR10502967.sra

2. sra格式转fastq格式

fastq-dump *.sra

二、数据质量过滤

01 检测数据质量

fastqc生成质控报告

fastqc *.fastq

multiqc将各个样本的质控报告整合为一个

multiqc *.zip

02 原始数据修剪

《原创10000+生信教程大神给你的RNA实战视频演练》中是双端测序，我以为自己分析的是单端测序，所以对代码稍作了修改，也不知道改得对不对。

运行如下代码，得到名为config的文件

mkdir $wkd/clean 
cd $wkd/clean 
ls /dat01/xietian/ncbi/public/sra/*.fastq >1
paste 1  > config

打开文件 qc.sh ，并且写入如下内容

bin_trim_galore=trim_galore
dir='/dat01/xietian/ncbi/public/sra/clean'
cat $1 |while read id
do
        arr=(${id})
        fq1=${arr[0]}
 $bin_trim_galore -q 25 --phred33 --length 36 --stringency 3 --fastqc -o $dir  $fq1
done

运行qc.sh

bash qc.sh config #config是传递进去的参数

三、比对到参考基因组

参考这篇文章：RNA-seq(5):序列比对：Hisat2

01 HISAT2官网下载index

cd /dat01/xietian/ncbi/public/sra
mkdir -p reference/index
cd /dat01/xietian/ncbi/public/sra/reference/index
wget -O hg38.tar.gz https://cloud.biohpc.swmed.edu/index.php/s/hg38/download
tar -zxvf *.tar.gz

02 开始比对：用hisat2，得到SAM文件

#进入目录
cd /dat01/xietian/ncbi/public/sra/aligned
#人的比对
#将以下内容写入aligned.sh文件：
for ((i=62;i<=67;i++));
do
hisat2 -t -p 8 -x /dat01/xietian/ncbi/public/sra/reference/index/hg38/genome -U /dat01/xietian/ncbi/public/sra/clean/SRR105029${i}_trimmed.fq -S SRR105029${i}.sam; 
done
#运行aligned.sh文件：
bash aligned.sh

只有3.25%比对上了吗，怎么肥四？？？

Time loading forward index: 00:00:03
Time loading reference: 00:00:00
Multiseed full-index search: 00:52:54
32061537 reads; of these:
  32061537 (100.00%) were unpaired; of these:
    31020728 (96.75%) aligned 0 times
    958065 (2.99%) aligned exactly 1 time
    82744 (0.26%) aligned >1 times
3.25% overall alignment rate
Time searching: 00:52:55
Overall time: 00:52:58

去网上搜索，指引我的又是Jimmy老师? 《做过1000遍RNA-seq的老司机告诉你如何翻车》（非常神奇啊，我踩过的坑，Jimmy老师）踩过

去NCBI上查了查，发现是PAIRED，误把双端测序当成了单端测序！

错误从这里开始（几乎是从头开始?)

2. sra格式转fastq格式

fastq-dump --split-files *.sra
#--split-files参数可以将其分解为两个fastq文件。
#如果不加该参数，则只有1个fastq文件（包含了两端测序的结果）。

转换得到的文件：

├── SRR10502962_1.fastq
├── SRR10502962_2.fastq
├── SRR10502963_1.fastq
├── SRR10502963_2.fastq
├── SRR10502964_1.fastq
├── SRR10502964_2.fastq
├── SRR10502965_1.fastq
├── SRR10502965_2.fastq
├── SRR10502966_1.fastq
├── SRR10502966_2.fastq
├── SRR10502967_1.fastq
└── SRR10502967_2.fastq

二、数据质量过滤

01 检测数据质量

fastqc生成质控报告，详细解读可以参考这篇文章?测序数据质量控制之FastQC

fastqc *.fastq

运行完成后得到这两种文件，打开html文件即可查看质控报告

*_fastqc.html
*_fastqc.zip

multiqc将各个样本的质控报告整合为一个

multiqc *.zip

02 原始数据修剪

运行如下代码，得到名为config的文件

mkdir $wkd/clean 
cd $wkd/clean 
ls /dat01/xietian/ncbi/public/sra/GSE140739/*_1.fastq >1
ls /dat01/xietian/ncbi/public/sra/GSE140739/*_2.fastq >2
paste 1 2  > config

打开文件 qc.sh ，并且写入如下内容

trim_galore，用于去除低质量和接头数据参数--fastqc：在数据过滤后再次质检

bin_trim_galore=trim_galore
dir='/dat01/xietian/ncbi/public/sra/GSE140739/clean'
cat $1 |while read id
do
        arr=(${id})
        fq1=${arr[0]}
        fq2=${arr[1]} 
 $bin_trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o $dir  $fq1 $fq2 
done

运行qc.sh

bash qc.sh config #config是传递进去的参数

感觉数据质量还不错，所以过滤这步这次被我忽略了。。。接着做

三、比对到参考基因组

由于测序仪机器读长的限制，在构建文库的过程中首先需要将DNA片段化，测序得到的序列只是基因组上的部分序列。为了确定测序reads在基因组上的位置，需要将reads比对回参考基因组上，这个步骤叫做mapping。目前mapping的工具有很多，这里用的是hisat2

01 HISAT2官网下载index

人类和小鼠的索引有现成的，HISAT2官网可以直接下载进行序列比对，这里下载的是人的index

cd /dat01/xietian/ncbi/public/sra
mkdir -p reference/index
cd /dat01/xietian/ncbi/public/sra/reference/index
wget -O hg38.tar.gz https://cloud.biohpc.swmed.edu/index.php/s/hg38/download
tar -zxvf *.tar.gz

02 开始比对：用hisat2，得到SAM文件

hisat2主要参数：

-x 参考基因组索引文件的前缀。 -1 双端测序结果的第一个文件。若有多组数据，使用逗号将文件分隔。Reads的长度可以不一致。 -2 双端测序结果的第二个文件。若有多组数据，使用逗号将文件分隔，并且文件顺序要和-1参数对应。Reads的长度可以不一致。 -U 单端数据文件。若有多组数据，使用逗号将文件分隔。可以和-1、-2参数同时使用。Reads的长度可以不一致。 –sra-acc 输入SRA登录号，比如SRR353653，SRR353654。多组数据之间使用逗号分隔。HISAT将自动下载并识别数据类型，进行比对。 -S 指定输出的SAM文件。

我的fastq文件在/dat01/xietian/ncbi/public/sra/clean

我的index在/dat01/xietian/ncbi/public/sra/reference/index/grch38

比对后得到的bam文件会存放在/dat01/xietian/ncbi/public/sra/GSE140739/aligned

#进入目录
cd /dat01/xietian/ncbi/public/sra/GSE140739/aligned
#人的比对
hisat2 -t -p 8 -x /dat01/xietian/ncbi/public/sra/reference/index/hg38/genome -1 /dat01/xietian/ncbi/public/sra/GSE140739/SRR10502962_1.fastq -2 /dat01/xietian/ncbi/public/sra/GSE140739/SRR10502962_2.fastq -S SRR10502962.sam

95.81%，终于成功了

Time loading forward index: 00:00:04
Time loading reference: 00:00:00
Multiseed full-index search: 00:29:13
32666485 reads; of these:
  32666485 (100.00%) were paired; of these:
    2361370 (7.23%) aligned concordantly 0 times
    28549069 (87.40%) aligned concordantly exactly 1 time
    1756046 (5.38%) aligned concordantly >1 times
    ----
    2361370 pairs aligned concordantly 0 times; of these:
      196235 (8.31%) aligned discordantly 1 time
    ----
    2165135 pairs aligned 0 times concordantly or discordantly; of these:
      4330270 mates make up the pairs; of these:
        2738848 (63.25%) aligned 0 times
        1444567 (33.36%) aligned exactly 1 time
        146855 (3.39%) aligned >1 times
95.81% overall alignment rate
Time searching: 00:29:14
Overall time: 00:29:18

03 SAM文件转BAM文件

samtools主要参数：

-view BAM-SAM/SAM-BAM转换和提取部分比对 -sort 比对排序 -merge 聚合多个排序比对 -index 索引排序比对 -faidx 建立FASTA索引，提取部分序列

for i in `seq 2 7`
do
#将sam文件转换成bam文件
    samtools view -S SRR1050296${i}.sam -b > SRR1050296${i}.bam
#对bam文件进行排序
#刚开始加了参数-o，运行后电脑终端一直不停跳乱码的东西，Xshell直接死机
    samtools sort SRR1050296${i}.bam SRR1050296${i}_sorted.bam
#对bam文件建立索引
    samtools index SRR1050296${i}_sorted.bam
done

运行完成可得到这三种文件：

SRR1050296${i}.bam
SRR1050296${i}_sorted.bam
SRR1050296${i}_sorted.bam.bai

此外还尝试过Jimmy老师的代码，但是报错了，后来又在生信技能树上找到了原因?samtools安装踩坑记录

04 加载IGV，可视化比对结果

载入参考序列，注释和BAM文件?学IGV必看的初级教程

四、基因表达水平分析

01 下载gtf基因组注释文件

用wget下载均失败了，最后用浏览器或者迅雷下载的，超级慢，还好只有三四十M

ftp://ftp.ensembl.org/pub/current_gtf/homo_sapiens/Homo_sapiens.GRCh38.99.gtf.gz
ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_29/gencode.v29.annotation.gtf.gz

02 使用featureCounts进行alignment-based的定量

featuresCounts软件用于统计基因/转录本上mapping的reads数，也就是用于raw count定量。该软件不仅支持基因/转录本的定量，也支持exon，gene bodies，genomic bins，chromsomal locations等区间的定量。详见?转录组定量可以用替换featureCounts代替HTSeq-count

featureCounts主要参数：

-a 输入GTF/GFF基因组注释文件 -p 这个参数是针对paired-end数据 -F 指定-a注释文件的格式，默认是GTF -g 从注释文件中提取Meta-features信息用于read count，默认是gene_id -t 跟-g一样的意思，其是默认将exon作为一个feature -o 输出文件 -T 多线程数

gtf="/dat01/xietian/ncbi/public/sra/reference/gtf/gencode/Homo_sapiens.GRCh38.99.gtf.gz" 
featureCounts -T 5 -p -t exon -g gene_id  -a $gtf -o  all.id.txt  *.bam
cat all.id.txt | cut -f1,7- > counts.txt #去除多余信息，保存表达矩阵为counts.txt

五、差异表达分析

从这里开始就是在R里面了，根据现成脚本稍微调试一下就OK，相对简单，就不具体写了

01 基因表达量的标准化方法及可视化

counts，RPKM，FPKM，TPM（?RNA-seq的counts值，RPM, RPKM, FPKM, TPM 的异同）

PCA图、热图等

02 差异表达分析及可视化

limma，voom，edgeR，DESeq2（?limma/voom，edgeR，DESeq2分析注意事项，差异分析表达矩阵与分组信息

差异基因的热图和火山图

03 三个软件包的差异分析结果比较及筛选

logFC 含义相关性图

实操花了整整两天，终于磕磕绊绊自己走了一遍，有种游戏玩通关的感觉?这次应该可以摆脱从入门到放弃的怪圈了，但代码的具体含义有些还要再理解，继续加油?

写在后面

活动链接：这120万我就不要了，送给500名优秀本科生，一切按照规则来，表明你足够优秀或者愿意拼命学习！就可以加入，生物信息学毕竟还是一个专业，大家作为生物学或者医学背景本科生，跨行学生物信息学技术，那么吃苦耐劳是基本素质！

还等什么，如果你想理解下面的图表，那就加入吧！

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2020-04-07，如有侵权请联系 cloudcommunity@tencent.com 删除

bash

bash 指令

本文分享自生信技能树微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

bash

bash 指令

登录后参与评论

0 条评论

热度