RNA-seq三代测序技术转录组上游分析

转录组RNA-SEQ上游分析2020

2020-09-05  本文已影响0人  superqun

安装配置conda

使用清华源下载sh脚本并安装

# 使用清华源下载sh脚本
wget -c  https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh

# 从官网下载最新版Miniconda3安装包，但速度较慢
wget -c https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh

下载完成后直接运行脚本文件bash Miniconda3-latest-Linux-x86_64.sh。需要输入yes然后等待安装完毕
最后安装好后，还不能马上使用conda，需要source一下bashrc

# 激活bashrc
source ~/.bashrc

==注意⚠️：==

设置镜像源

# 下面这四行配置清华大学的bioconda的channel地址，国内用户推荐
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
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 --add channels r 
conda config --add channels conda-forge 
conda config --add channels bioconda

在设置后镜像或者设置不自动进入base后，会在.condarc文件中自动生成config信息。如下：

$ cat .condarc 

auto_activate_base: false
channels:
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
  - defaults

conda环境创建

创建一个python2的环境管理：

conda create -y -n rna_seq python=3

# -y        自动确认
# -n        新环境名字
# python=3  新环境中python=3

激活和退出环境

conda activate <conda_name>     #激活某环境
conda decativate <conda>        #取消激活某环境

conda安装软件

在软件环境中使用命令安装软件

conda install -y sra-tools      #安装sra-tool软件，可以通过空格安装多个软件
conda install -y sra-tools fastqc trim-galore hisat2 subread multiqc samtools salmon fastp

conda软件安装位置和普通软件安装位置不一样，通过which <softname>来查看conda安装的软件位置

质量评估 @ ==fastQC==

fastq格式

FastQ格式描述：https://mp.weixin.qq.com/s/8g-oUjiEhV4cGMJNuhmISQ
FastQ格式wiki：https://en.wikipedia.org/wiki/FASTQ_format
FastQ格式文献：https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2847217/

概念
FastQ格式是序列格式中常见的一种，它存储了生物序列以及相应的质量评价，其序列以及质量信息都是使用一个ASCII字符标示，最初由Sanger开发，目的是将FASTA序列与质量数据放到一起，目前已经成为高通量测序结果的事实标准。

格式说明
FASTQ文件中每个序列通常有四行：

  • 1.第一行：必须以“@”开头，后面跟着唯一的序列ID标识符，然后跟着可选的序列描述内容，标识符与描述内容用空格分开；

FsatQC软件

FastQC质量评估软件官网：http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
==注意⚠️==

常用参数：

# 常用参数
fastqc -o <out.dir> -t <thred_num> -f <input_format>  <input_file_1> <input_file_2> ...

# -o    设置输出目录
# -t    设置线程数
# -f    设置输入文件格式

批处理

# bash中
a=`ls *.fq`
fastqc -o ./fastqc_raw -t 10 $a

# zsh中
b=`ls -C *.fq`
fastqc -o ./fastqc_raw -t 10 `echo $b`

multiqc综合所有qc

使用multiqc来把所有的质量评估放在一起观察

multqc -o <out.path> *.fastqc.zip

结果解读

单一碱基占比

image

unique reads 的占比，

单碱基测序质量在150bp上的分布情况

image

也有这样的

image

当前RNA-seq测序技术，测序错误率分布存在以下两个特征。

  1. 测序错误率随着测序序列(Sequenced Reads)长度的增加而升高。这是由测序过程中化学试剂的消耗导致的，为Illumina高通量测序平台所具有的特征。
  2. 前6个碱基具有较高的测序错误率，此长度恰好为RNA-seq建库过程中反转录所需的随机引物长度。前6个碱基测序错误率较高是因为随机引物和RNA模版的不完全结合(Jiang et al.)。

测序质量的分布图

image

和单碱基测序质量在150bp上的分布情况不同，这个是单个样本中，碱基质量的分布情况。绝大部份集中在Q30以上，效果良好（类似于密度分布图）

GC含量测定

1.整体GC含量测定，主要看是否有双峰，如果有双峰可能有其他物种掺入。（GC含量在物种间存在一定特异性）


image

2.单碱基GC(TAN)在150bp上的分布情况。理想情况是四条线在25%轻微波动，但是如所见，前几个bp非常不稳定。这是由于反转录过程中所使用的6bp随机引物，会引起前几位碱基在核苷酸组成上有一定偏好性，产生正常波动，随后则趋于稳定。

image
image

下图除了加入了ATCG等碱基，还加入了不确定碱基N碱基的含量。后期会去除这些含N碱基的reads。

接头adapter统计

image

统计了接头出现的累积频率。一般累积频率不超过5%即可

过滤reads @ ==trim_galore==

过滤标准

一般的过滤标准：

trim_galore

trim_galore使用手册：https://github.com/FelixKrueger/TrimGalore/blob/master/Docs/Trim_Galore_User_Guide.md
阅读延伸：https://www.jianshu.com/p/7a3de6b8e503

Trim Galore是对FastQC和Cutadapt的包装。适用于所有高通量测序，包括RRBS(Reduced Representation Bisulfite-Seq ), Illumina、Nextera 和smallRNA测序平台的双端和单端数据。主要功能包括两步：
第一步首先去除低质量碱基，然后去除3' 末端的adapter, 如果没有指定具体的adapter，程序会自动检测前1million的序列，然后对比前12-13bp的序列是否符合以下类型的adapter:

==注意⚠️==

# 常用参数
trim_galore

-j                  #使用线程数,注意假设已使用Python 3并安装了Pigz，那么内核设置为4，实际使用内核是15，因此，最高设为4
-q                  #切除质量得分低于设定值的序列
--phred33           #得分标准（默认）
-a                  #输入adapter序列，也可以不输入，软件自动搜索可能性最高的平台对应的adapter，也可以直接设置参数输入这三个平台分别是--illumina --nextera --small_ma
--illumina          #设置测序平台是ilumina，用来设定接头序列用
--length            #设置小于此bp的reads会被删除，注意，在pe150下，不要设置太高，可以50或36（默认20）
--max_length        # 设置长度大于此值被丢弃
-s/--stringency     #设定可以容忍的前后adapter的重叠碱基数，默认是1，可以放宽以为后一个adapter几乎测不到
--paired            #上端reads，如果一个被剔除，则另一个也被剔除
-o                  #设定已经存在的输出目录
--fastqc            #当分析结束后，使用默认选项对结果文件进行fastqc分析
--gzip              #输出结果使用gzip压缩,gzip压缩将大大减慢多核进程，因此几乎不值得
--max_n             #设置如果超过设定n就删除。可以根据length选项，
--basename          #设置输出文件的基本名称，尽在使用1个文件或2个文件参数时候可用

例子

trim_galore --phred33 -q 20 --length 36 --stringency 3 --fastqc --paired --max_n 3 -j 4 -o ./trim_galore 1.fastq.gz 2.fastq.gz

批处理

ls *1.fq >1.txt
ls *2.fq >2.txt
paste 1.txt 2.txt >config

# cat trim_galore.sh
cat config |while read id
do
    arr=(${id})
    fq1=${arr[0]}
    fq2=${arr[1]}
    echo "trim_galore work_ID：${fq1} started at $(date)"
    echo "trim_galore -q 20 --phred33 --stringency 3 --length 36 --paired --fastqc  --max_n 3 -j 4 ${fq1} ${fq2}"
    echo "trim_galore work_ID：${fq1} finished at $(date)"
    echo -e "\n"
done

# 进行上述脚本
nohup sh trim_galore.sh >trim_galore.log &

批处理参考：https://www.jianshu.com/p/5a77b4156a5d
并行处理要求：https://www.jianshu.com/p/74a328f4d23a

比对 @ ==hisat2==

hisat用户手册：http://ccb.jhu.edu/software/hisat2/manual.shtml
==简书：有关hisat深入介绍，推荐==：https://www.jianshu.com/p/ce3f4afb9b60

hisat的优势：

参考基因组，gtf文件选择，不同基因组版本差异

==注意⚠️：==

参考基因组和gtf文件可以从esemble下载。资源汇总地址：ftp://ftp.ensembl.org/pub/ 。可以选择最新的 101版本（截止2020/08/25）。

image
在101版本中选gtf来下载最新的gtf注释文件，选择fasta来选择最新的基因组文件。
基因组文件依次选择release-101 > fasta > homo_sapiens > dna 然后在诸多文件中选择*GRCH38.dna.primary_assembly*
image

不同基因组版本差异：

参考&资料：

构建索引

人的基因组索引可以下，但是可能会有基因组差异问题，所以，除非非常清楚和索引的基因组匹配，否则建议自己构建索引。当然，构建索引是一个非常漫长的过程。

hisat的索引策略

hisat的比对策略

image

索引的构建

注意，如果要用到gtf中exon和split site信息，要先用hisat2自带的py软件生成，然后添加到参数--exon 和参数--ss后

hisat2-build [options] <reference_in> <ht2_index_base>
# 常用[options]：
# -p                线程数
# 参数：
# <reference_in>    参考基因组fasta文件 list，如果为list，使用逗号分开
# <ht2_base>        索引文件的前缀名
# --ss              > 与 --exon 一起使用，提供拼接位点信息。只支持hisat2自己格式
#                   ^（需要用hisat2_extract_splice_sites.py 与 gtf来生成）
# --exon            指定外显子列表（只支持hisat2自己格式，需要通过hisat2_extract_exons.py 与gtf文件来生成）从GTF提取剪接位点。


# 例子
nohup \ 
hisat2-build -p 6 Homo_sapiens.GRCh38.dna.primary_assembly.fa homo_GRCh38.dna.primary \
&

索引构建完毕如下：

image

hisat2比对

# 常用参数
hisat2 [option] -1 <m1> -2 <m2> -x <index_base> -S <out.sam>

-x <path/to/base>             #索引的前缀,也可以添加路径，例如 ～/genome/homo_GRCh38.dna.primary
-1              #双端测序的第一个文件，若有多组，则用逗号隔开，名字与2要匹配，如-1 flyA_1.fq,flyB_1.fq
-2              #类比上，名字与1要匹配，如-2 flyA_2.fq,flyB_2.fq
-U              #单端数据文件
-S              #保存到的sam文件，默认输出到标准输出

# 以下为参数
-p              #线程数
--dta           # > 在下游使用stringtie组装的时候量身定做的参数。使用此选项，
                # ^ HISAT2需要更长的锚长度才能从头发现拼接位点，这样可以减少与短锚的对齐，
                # ^ 这有助于转录汇编程序显着提高计算和内存使用率。
                # ^ 当然下游不使用stringtie也可以使用此参数减少短锚定read
--dta-cufflinks #下游使用cufflinks则需要添加此参数
-q              #指定读取的测序文件是fastq文件（含有质量信息）
-f              #指定读取的测序文件是fa文件
-t              #打印加载索引文件和对齐读数所需的时钟时间
--phred33       #质量表示方法，如果使用fq文件，建议选择，同理有小众的--phred64
--min-intronlen #设置最小内含子的长度，默认值20
--max-intronlen #设置最大内含子长度，默认500000
--known-splicesite-infile <path>    #提供已知的剪接位点列表，HISAT2利用这些列表比对,可以用hisat2_extract_splice_sites.py和gtf生成信息
--novel-splicesite-outfile <path>   #在结果中生成一个剪接位点的列表
--novel-splicesite-infile <path>    # > 可以使用novel-splicesite-outfile所生成的剪接列表作为
                                    # ^ 新剪接列表，应该可以自己定义。为特定基因。

例子

# 例子 
hisat2 -p 10 -q --known-splicesite-infile <splicesites.txt> -x <homo_GRCh38_dna_primary> -1 <S_CTRL_1_val_1.fq.gz>,<S_OE_1_val_1.fq.gz> -2 <S_CTRL_2_val_2.fq.gz>,<S_OE_2_val_2.fq.gz> -S <./sam/sample.hisat.sam>

批处理

ls *1.fq.gz|cut -d"_" -f1,2>id.txt      #整理列举id，在trim_galore目录下，用已经过滤好的数据

# cat hista_align.sh
cat ~/lab/ZhangJinRui/analysis/trim_galore/id.txt | while read id                                
do
    echo "hista2 :work_ID ${id} start at $(date)"
    echo "hisat2 -p 10 -q --known-splicesite-infile ~/lab/ZhangJinRui/analysis/geom/splicesites.txt -x homo_GRCh38_dna_primary -1 ~/lab/ZhangJinRui/analysis/trim_galore/${id}_1_val_1.fq.gz -2 ~/lab/ZhangJinRui/analysis/trim_galore/${id}_2_val_2.fq.gz -S ~/lab/ZhangJinRui/analysis/sam/${id}.hisat.sam"
    echo "hisat2 :work_ID ${id} finished at $(date)"
done

# 运行上述脚本 在索引目录下？不用也可以，索引前可以用路径指定 -x 
nohup sh hista_align.sh >hisat_align.log &

PS：
使用hisat2_extract_splice_sites.py来生成gtf中的剪接位点，辅助比对。使用方法：

hisat2_extract_splice_sites.py genes.gtf > splicesites.txt

sam格式转换 @ ==samtools==

sam格式：The Sequence Alignment / Map format，存储了序列比对情况的文件，是一个文本文件，都以tab分列。体积较大。分为头部区和主体区。
官方文档：http://samtools.github.io/hts-specs/SAMv1.pdf

头部区

image

其他具体sam格式信息：https://blog.csdn.net/genome_denovo/article/details/78712972

PS:诺和分析数据的代码：

/TJPROJ6/RNA_T/software/miniconda2/envs/hisat2-2.0.5/bin/hisat2-align-s --wrapper basic-0 -x /TJPROJ6/GB_TR/reference_data/new_pip/Animal/Homo_sapiens/Homo_sapiens_Ensemble_94/Homo_sapiens_Ensemble_94 -p 4 --dta -t --phred33 --passthrough -1 /tmp/22775.inpipe1 -2 /tmp/22775.inpipe2

排序和转BAM

samtools官方文档：http://www.htslib.org/doc/samtools.html
==注意⚠️：==

排序转bam

详细使用方法可以使用man samtools sort查看

samtools sort [option] <input.sam>
            -@                      #线程数
            -o <path/to/out.bam>    #输出的sort后文件
            -T <path/to/temp>       #临时文件
            -O <sam|bam..>          #输出文件格式，默认根据-o来确定
            -n                      #按照QNAME字段排序而不是染色体坐标
            
# 例子
samtools sort -@ 10 -o ./sort_bam/sample_sort.bam sample.sam

批处理

ls *.sam|cut -d"." -f1,2 >sam.id.txt        #把诸如 S_CTRL.hisat.sam S_OE.hisat.sam 取前缀名然后分行存储在文本内

# cat sam_sort.sh 
cat sam.id.txt|while read id  
do
    echo -e "*******WORK_ID ${id} START AT : \n$(date)"
    samtools sort -@ 10 -o ./sort_bam/${id}.sort.bam ${id}.sam
    echo -e "*******WORK_ID ${id} FINISHED AT : \n$(date)"
done

# 运行脚本 在当前sam文件目录下，并且在目录下新建了 sort_bam文件
nohup sh sam_sort.sh >sam_sort.log &

索引bam文件

samtools index sample.sort.bam  #生成文件默认是原文件名 + .bai 也可以指定在原文件后面直接指定

批处理

ls *.bam > bam.list.txt

#cat bam_index.sh
cat ba.list.txt|while read id
do
    echo -e "*******WORK_ID ${id} START AT : \n$(date)"
    samtools index ${id}
    echo -e "*******WORK_ID ${id} FINISHED AT : \n$(date)"
done 

#运行
nohup sh bam_index.sh >bam_index.log &

==注意⚠️：==
说明网站上指出，samtools index 有一个 -@ 指定线程的参数。但是自己实际操作报错。

计数 @ ==featureCounts==

英文简明说明书：http://bioinf.wehi.edu.au/featureCounts/
官方网站：http://subread.sourceforge.net/
英文用户手册下载：http://bioinf.wehi.edu.au/subread-package/SubreadUsersGuide.pdf （featurecounts说明在P31）
简书详细介绍（推荐）：https://www.jianshu.com/p/9cc4e8657d62

featureCounts理解的核心是feature。我们感兴趣的基因组特征（feature）可能是外显子、基因、启动子区域、gene body或者是基因组间隙，不单单是gene_id。

在计数中RNA-seq计数可能更加复杂，：

输入的GTF的作用：
GFF/GTF/SAF主要提供feature identifier(如geneID), chromosomename, start position, end position and strand 这五列信息

gtf信息解读：
https://www.jianshu.com/p/7cec5a02768a

feature和meta-feature：

read和fragment：
如果是双端测序，这一对read定义了一个DNA/RNA片段的两端，这种情况下，featurecounts会计算片段数（fragment）而不是read数。即一个fragment由两个raed确定。（？fragment可能包含有read中没有的bp信息，但是read负责确定fragment的定位信息）

多重overlap的选择：

==注意⚠️：==

主要参数

featureCounts [options] <input.file>

     <input.file>       #输入文件，sam/bam 支持多个文件输入
     -a <gtf>           #参考gtf注释文件，支持gz压缩
     -F <参考文件格式>    #默认GTF
     -T <int>           #线程数 1-32
     -J                 #对可变剪接计数
     -G <str>           #有-J设置时，用来提供参考基因组辅助寻找可变剪接
     -M                 #如果设置了-M则多重比对会被统计
     -o <out.file>      #输出文件的名字，输出文件的内容是read的统计数目
     -O                 #允许多重比对，当一个read比对到多个feature或多个metafeature的时候，这条read会被统计多次
     -p                 #声明测序类型是paired-end，此时，后续统计fragment而不是read
     -B                 #在-p设置时，只有两端都比对上的fragment才会被统计
     -C                 #在-p设置时，-C声明比对到不同染色体的fragment不计数
     -d <int>           #最短的fragment ，默认50
     -D <int>           #最长的fragment，默认600
     -f                 #设置后会统计feature水平数据，如exon-level，否则会统计meta-feature层面数据，如gene-level，（？应该是和-g冲突,一般与-t连用?）
     -g <str>           #参考的gtf提供的时候，我们需要提供一个id identifier 来将feature水平的统计汇总为meta-feature水平的统计，默认为gene_id，可以选择transcript_id、gene_name、transcript_name等
     -t <str>           #设置feature-type，-t指定的必须是gtf中有的feature，同时read只有落到这些feature上才会被统计到，默认是“exon”（？应该和-g冲突）

例子：

featureCounts -T 10 -a $gtf -o all.id.txt -p -t exon -g gene_id   *.sorted.bam

featureCounts -t exon -g gene_id -a annotation.gtf -o counts.txt library1.bam library2.bam library3.bam

featureCounts -p -t exon -g gene_id -a annotation.gtf -o counts.txt mapping_results_PE.bam

featureCounts -p -C -t exon -g gene_id -a annotation.gtf -o counts.txt mapping_results_PE.bam

featureCounts -p -B -t exon -g gene_id -a annotation.gtf -o counts.txt mapping_results_PE.bam

#输出gene_id
nohup featureCounts -T 10 -a ~/lab/ZhangJinRui/analysis/geom/Homo_sapiens.GRCh38.101.gtf -o counts.txt -p -t exon -g gene_id ~/lab/ZhangJinRui/analysis/sam/sort_bam/*.sort.bam &

#输出gene_name
nohup featureCounts -T 10 -a ~/lab/ZhangJinRui/analysis/geom/Homo_sapiens.GRCh38.101.gtf -o gene_name_counts.txt -p -t exon -C -B -g gene_name ~/lab/ZhangJinRui/analysis/sam/sort_bam/*.sort.bam >featurecount_gene_name.log &

结果处理

一般输出结果大于7列，其中第一列是由-g决定的meta-feature，第七列开始是对应bam文件的count数，第七列后面多样本的couts数，如果要引用多个数，例如，共20个样本，要选从第七列选到26列，可以使用

cat counts.txt|head|cut -f1,`seq -s"," 7 26`

我想建立并管理一个高质量的生信&统计相关的微信讨论群，如果你想参与讨论，可以添加微信：veryqun 。我会拉你进群，当然有问题也可以微信咨询我。









上一篇下一篇

猜你喜欢

热点阅读