HiCUP的应用以至结果解读,Seq数据深入分析

9778818威尼斯官网 1

SRA—>FASTQ—>BAM—>COUNTS 这几个步骤而已,中间穿插一些质控的手段,每个步骤选择好合适的软件即可。可以参考:一个植物转录组项目的实战 http://www.bio-info-trainee.com/2809.html

使用的数据,双端测序数据clean data 1 & clean data 2**

map与alignment

新一代测序数据的处理过程中,map/alignment总是同时出现,那么二者的区别是什么呢?大致上可以这样理解:map这个动作表示的是:把短片段对应到参考基因组的某个位置上。这个动作只回答“短片段来自长片段的什么位置?”这个问题;而align这个动作早在二代测序之前很久就有了,这个动作表示两条(或多条)序列(不管是核酸序列还是蛋白序列)的比较,回答的问题是“这两条序列之间有什么异同?”。显然在处理重测序数据的过程中,要先map,再align,即先把read定位到参考基因组的某个位置上,再比较read和参考基因组序列的异同。另外还可以认为:凡是要align,那一定是要先map的;但是map过后,却不一定要align(比如我只关注CNV或表达量等信息,而不关心具体序列的时候,就只需要map确定下read来自参考序列的哪个位置,而不需要align获取具体的序列间异同信息)。可见在SAM文件中,map和align这两个动作都要有,这可能就是这两个动作总是同时出现的原因。理解了这两个名字的含义和区别,就可以理解为什么是“mapping” position,因为确定position信息只需要mapping。map和align两个词,似乎都可以用中文“比对”一词来表示。下文中如果没有必要区分二者,将全都译作“比对”。

通过前面两篇文章对MutMap原理和分析流程的介绍,我们对MutMap的分析已经有了比较深入的了解,是时候来点儿数据实战一下了!

HiCUP需依赖Perl,Bowtie/Bowtie2,R以及Samtools

这是RNA-Seq 上游分析的大致流程,比对 定量。当然实验目的若只需要定量已知基因,也可以选择free-alignment 的流程工具如kallisto/Salmon/Sailfish,其优点是可用于RNA-seq的基因表达的快速定量,但是对于小RNA和表达量低的基因分析效果并不好(2018年刚发表的一篇文章对free-alignment 的工具进行了质量评估,doi: https://doi.org/10.1101/246967)。基于比对的流程,比对工具也有很多选择,如Hisat,STAR,Tophat(hista可以替代tophat),bowtie等, 还有据说速度超快的Subread。根据比对工具的选择,定量软件也可以配套选择,如
STAR/bowtie—RSEM;
Hisat —featureCounts/HTseq;
Subread—featureCounts。
由于后续要鉴定新的转录本,因此需要对转录本组装,组装可以选择cufflinks,或者stringtie。文章是用的Tophat cufflinks。我的后续分析打算用Hisat2比对,stringtie组装,featureCounts定量,同时打算尝试下Subread featureCounts的流程。
Hisat2 stringtie featureCounts;
Subread featureCounts

1、bowtie建库

bowtie2-build genome.fa genome_index (*.fa代表 要建库的文件名字,genome_index代表建库后的名字)


bwa

MutMap的发明者已经搭建好了一套Pepline,并写好了MutMap pipeline framework供大家免费使用。我们只需要将数据按照要求进行放置和命名,并通过非常简单的几个命令,将任务提交给服务器,经过三五天的分析(以水稻为例),就可以拿到结果了。

在官网


2、将read比对到基因组中

bowtie    -q -S -t  -p 20   -m 1   --best --strata  genome_index input.fastq  out.sam

-q 指输入的测序数据是Fastq格式

-S 输出的结果文件格式

-p 使用的线程数

-m 只保留有一个alignment的read

命令算法:

bwa <command> [options]

bwa index 建立索引
bwa mem BWA_MEM algorithm
bwa bwasw BWA_SW for long queries
bwa aln gapped/ungapped alignment
bwa samse generate alignment(single ended)
bwa sampe generate alignment(paired ended)

可通过bwa index/bwa mem/bwa aln查看相应命令的帮助


wget

SRA—>FASTQ

sra格式的数据需要先用fastq-dump转换, --split-3 表示双端测序,--gzip将生成的fastq文件压缩。

cd /pnas/fangxd_group/renyx/macaque/01rawdata
for ((i=230;i<=237;i  )) ;do /software/biosoft/software/sratoolkit.2.8.2-1-centos_linux64/bin/fastq-dump --split-3 --gzip SRR4042$i.sra;done

3、查看比对结果

samtools view -S -b sample.sam | samtools flagstat -


比对流程

-建立一个已知参考基因组的索引

索引的建立是对参考基因组进行重新格式化,每个比对软件的索引不同。基因组越大,构建索引的时间越长。很多比对工具在其官方网站上有提供已经预建的索引文件,可直接下载使用。

-FASTA/FASTQ文件中的reads比对到索引上

1.分析环境的搭建

HiCUP的应用以至结果解读,Seq数据深入分析。首先,我们需要有一台Linux服务器。由于分析过程中需要进行大量的运算,个人PC的配置根本吃不消。所以我们需要一台服务器,当然配置越高越好。在服务器上安装Linux系统,我们使用的是Redhat,当然也可以选择免费的Ubuntu或CentOS。

然后,在服务器安装如下软件:

Perl (v5.8.8) 5.22.2

R (version 2.15.0)3.0.1

BWA (version 0.5.9-r16)本版

SAMtools (0.1.8 or before) 本版

FASTX-Toolkit 使用conda安装0.0.14版

(http://hannonlab.cshl.edu/fastx_toolkit/download.html)

(括号中的为推荐版本,后面的黑体为我自己使用的版本)

关于FASTX-Toolkit版本的选择:

一定要安装0.0.13及以上的版本,因为从0.0.13开始FASTX-Toolkit才开始支持Illumina 1.8 Phred 33的质量打分规则,否则会报错。就这一个报错,花费了我一天的时间,直到使用conda安装了0.0.14版本才解决。(2018.01.24记)。

# 下载到压缩包hicup_v0.6.1.tar.gz

质控

mkdir QC
echo "fastqc started at $(date)"     #输出运行时间
fastqc -o QC *gz
multiqc *fastqc.zip --ignore *.html
echo "fastqc finished at $(date)"

质控检测采用fastqc和multiqc联合使用, multiqc有利于多个样本同时展示质量检测结果。FastQC的检测报告左侧会给出测序质量的一个summary,通常并不是所有的检测项都会是绿色通过,会有一些警告和fail的项目。主要可以从Per base sequence quality 看一下测序碱基质量,Per sequence GC content 看一下GC含量,如果实际的GC含量(红线)出现双峰,且导致后期的序列比对很低时,需要引起注意,有可能是存在样品污染;再者就是看一下Adapter是否去除及去除的是否干净。这里的Adapter虽然显示通过,但是去除的并不干净,所以后续还会进一步去除Adapter。

9778818威尼斯官网 2

Summary QC

9778818威尼斯官网 3

Adapter Content

但为了使用后面的cufflink套件,一般使用tophat做序列比对

示例

参照原文用ebola病毒测序数据和1976年Mayinga参考基因组。
mkdir -p ~/refs/ebola #建立文件夹用于存放参考基因组文件及将建立的索引存放于此
下载ebola病毒参考基因组文件,命名为1976.fa
efetch -db=nuccore -format=fasta -id=AF086833 > ~/refs/ebola/1976.fa
建立索引
bwa index ~/refs/ebola/1976.fa
还好前段时间看了下shell编程,才看懂文中的变量的用法,基础差的学生就得想办法恶补。

REF=~/refs/ebola/1976.fa #将后面的路径赋值给变量REF
bwa index $REF #shell中用$加变量名用来引用变量值
ls ~/refs/ebola #查看建立索引后目录内文件的变化

1976.fa 1976.fa.amb 1976.fa.ann 1976.fa.bwt 1976.fa.pac 1976.fa.sa

获取SRR1972739的ebola病毒的测序数据

esearch -db sra -query PRJNA257197  | efetch -format runinfo > runinfo.csv
fastq-dump -X 10000 --split-files SRR2571197#通过fastq-dump命令将sra数据还原为fastq文件,双端测序文件加上参数--split-files

命令结束后会在当前目录下产生两个测序数据
SRR1972739_1.fastq和SRR1972739_2.fastq

R1=SRR1972739_1.fastq #定义定量R1
R2=SRR1972739_2.fastq #定义变量R2

执行比对
bwa mem $REF $R1 $R2 > bwa.sam

bwa mem

9778818威尼斯官网 4

image.png

2.MutMap pipeline framework的下载和安装

先下载MutMap pepline程序,解压到Linux系统的工作目录下(此处以myhome为例)。

#切换到安装目录

$ cd /myhome

#将framework压缩文件拷贝到相应文件夹

$ cp MutMap_framework1.4.4.tar.gz

#解压

$ tar zxvf MutMap_framework1.4.4.tar.gz

#改名

$ mv MutMap_framework1.4.4 MutMap_test

tar -zxvf hicup_v0.6.1.tar.gz

去除接头和低质量值

Trimmomatic 支持多线程,处理数据速度快,主要用来去除 Illumina 平台的 fastq 序列中的接头,并根据碱基质量值对 fastq 进行修剪。软件有两种过滤模式,分别对应 SE 和 PE 测序数据,同时支持 gzip 和 bzip2 压缩文件。另外也支持 phred-33 和 phred-64 格式互相转化。

2、使用Tophat做序列比对

依赖软件有bowtie、samtools,不同版本的Tophat对bowtie的要求不同

tophat -p 8 -G annotation.gtf  -o tophat_output bowtie_index pair_end-1 pair_end-2

-p 使用的线程数

-G 使用的注释文件 

-o 输出文件夹

注意 1、bowtie_index的名字尽量和基因组的名字(不包括后面的.fa)一样

         2、此处使用的基因组和注释文件必去相互对应,即基因组各个染色体的名字必须与注释文件的第一列对应

bowtie

Bowtie 2 is an ultrafast and memory-efficient tool for aligning sequencing reads to long reference sequences. It is particularly good at aligning reads of about 50 up to 100s or 1,000s of characters, and particularly good at aligning to relatively long (e.g. mammalian) genomes. Bowtie 2 indexes the genome with an FM Index to keep its memory footprint small: for the human genome, its memory footprint is typically around 3.2 GB. Bowtie 2 supports gapped, local, and paired-end alignment modes.

算法可参考天天向上的博客

3.将测序数据链接到相应的文件夹中

将数据上传到服务器中专门的目录下,通过链接的方式,将原始文件链接到相应的操作目录下。

# HiCUP是Perl写的解压即可使用,无需编译

Trimmomatic 过滤步骤

Trimmomatic 过滤数据的步骤与命令行中过滤参数的顺序有关,通常的过滤步骤如下:
ILLUMINACLIP: 过滤 reads 中的 Illumina 测序接头和引物序列,并决定是否去除反向互补的 R1/R2 中的 R2, 该引物序列可以在Trimmomatic软件的安装目录下找到,双端通常选择TruSeq3-PE-2。
SLIDINGWINDOW: 从 reads 的 5’ 端开始,进行滑窗质量过滤,切掉碱基质量平均值低于阈值的滑窗。
MAXINFO: 一个自动调整的过滤选项,在保证 reads 长度的情况下尽量降低测序错误率,最大化 reads 的使用价值。
LEADING: 从 reads 的开头切除质量值低于阈值的碱基。
TRAILING: 从 reads 的末尾开始切除质量值低于阈值的碱基。
CROP: 从 reads 的末尾切掉部分碱基使得 reads 达到指定长度。
HEADCROP: 从 reads 的开头切掉指定数量的碱基。
MINLEN: 如果经过剪切后 reads 的长度低于阈值则丢弃这条 reads。
AVGQUAL: 如果 reads 的平均碱基质量值低于阈值则丢弃这条 reads。
TOPHRED33: 将 reads 的碱基质量值体系转为 phred-33。
TOPHRED64: 将 reads 的碱基质量值体系转为 phred-64。
参考:NGS 数据过滤之 Trimmomatic 详细说明

echo "trimmomatic cut adapters started at $(date)"
for i in {230..237}
do
java -jar /software/biosoft/software/Trimmomatic-0.36/trimmomatic-0.36.jar PE -threads 4 SRR4042$i_1.fastq.gz SRR4042$i_2.fastq.gz 
SRR4042$i_paired_clean_R1.fastq.gz 
SRR4042$i_unpair_clean_R1.fastq.gz 
SRR4042$i_paired_clean_R2.fastq.gz 
SRR4042$i_unpair_clean_R2.fastq.gz 
ILLUMINACLIP:/software/biosoft/software/Trimmomatic-0.36/adapters/TruSeq3-PE-2.fa:2:30:10:1:true 
LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:50 TOPHRED33
done
mv *_paired_clean* ../02clean_data/
echo "trimmomatic cut adapters finished at $(date)"

<ATTENTION>

         Mapping的最后一步是去除map到基因组中多于一处的序列,如果出现好几个序列都map在完全相同的一个区段,那么就应该只保留一个这样的序列,所以,只保留匹配最高的那一个。而且这样的序列占很大一部分,这步也很简单,samtools里的rmdup可以轻松解决:

samtools rmdup -s input.bam output.bam

工作流程

-建立索引
-将FASTA/FASTQ格式的reads与索引进行比对

bowtie2 -h

Usage:
bowtie2 [options]* -x <bt2-idx> {-1 <m1> -2 <m2> | -U <r>} [-S <sam>]

<bt2-idx> Index filename prefix (minus trailing .X.bt2).
NOTE: Bowtie 1 and Bowtie 2 indexes are not compatible.
<m1> Files with #1 mates, paired with files in <m2>.
Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
<m2> Files with #2 mates, paired with files in <m1>.
Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
<r> Files with unpaired reads.
Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
<sam> File for SAM output (default: stdout)

<m1>, <m2>, <r> can be comma-separated lists (no whitespace) and can be
HiCUP的应用以至结果解读,Seq数据深入分析。specified many times. E.g. '-U file1.fq,file2.fq -U file3.fq'.

建立索引文件
bowtie2-build -h

Usage: bowtie2-build [options]* <reference_in> <bt2_index_base>
reference_in comma-separated list of files with ref sequences
bt2_index_base write bt2 data to files with this dir/basename
*** Bowtie 2 indexes work only with v2 (not v1). Likewise for v1 indexes. ***
Options:
-f reference files are Fasta (default)
-c reference sequences given on cmd line (as
<reference_in>)
--large-index force generated index to be 'large', even if ref
has fewer than 4 billion nucleotides
-a/--noauto disable automatic -p/--bmax/--dcv memory-fitting
-p/--packed use packed strings internally; slower, less memory
--bmax <int> max bucket sz for blockwise suffix-array builder
--bmaxdivn <int> max bucket sz as divisor of ref len (default: 4)
--dcv <int> diff-cover period for blockwise (default: 1024)
--nodc disable diff-cover (algorithm becomes quadratic)
-r/--noref don't build .3/.4 index files
-3/--justref just build .3/.4 index files
-o/--offrate <int> SA is sampled every 2^<int> BWT chars (default: 5)
-t/--ftabchars <int> # of chars consumed in initial lookup (default: 10)
--threads <int> # of threads
--seed <int> seed for random number generator
-q/--quiet verbose output (for debugging)
-h/--help print detailed description of tool and its options
--usage print this usage message
--version print version information and quit

bowtie2-build $REF $REF #建立索引

第一个路径表示参考文件,第二个表示索引的前缀名称,可以是anything。

ls $REF.*
运行比对
bowtie2 -x $REF -1 $R1 -2 $R2 > bowtie.sam

3.1 通过MD5值验证文件的完整性

HiCUP中的比对工具是Bowtie/Bowtie2,建议采用最新版的Bowtie2。和HiCUP一样,Bowtie2有Linux二进制版,解压即可使用,无需编译

质控前后的结果

9778818威尼斯官网 5

General Statistics


3、使用cufflink套件计算表达量和差异表达基因

Cufflinks是一套拼接转录本,计算表达量,计算差异表达的工具。尽可能拼接处最优可能的转录本的结构,并且估计它的表达量。

cufflinks -p 10 -u -g/G annotation.gtf -o Cufflink1_R1_output tophat_output/accepted_hits.bam

-p 使用的线程数

-u 告诉Cufflinks用更准确的方法去处理贴到多个位点上的reads,如果没有-u,Cufflinks只会把这些reads简单地平均分配。比如一个read贴到了10个位置,那么每个位置分得十分之一。加-u后会先进行平均分配,然后按照这10个位置各自的表达量,计算read被分配到每个位置的概率。

-G是提供一个注释文件,并且告诉Cufflinks不要去拼接新的转录本,只能用注释文件里提供的转录本。

-g 也是提供一个注释文件,但是Cufflinks会在这些已知转录本的指导下,拼接新的转录本

-o 输出的文件夹

最后加的bam文件是Tophat运行结束产生的结果

比对软件的比较与选择

(1)Windows中生成MD5值

参考链接:https://jingyan.baidu.com/article/f7ff0bfc18ff302e26bb1382.html

  1. HiCUP的比对策略

Hisat2 比对 featureCounts定量

Hisat2 比对,首先用hisat2-build 构建index,然后比对,输出sam格式,再用samtools将sam转为bam格式,并排序构建索引。
featureCounts 目前已经整合在subread,subread是一个综合的软件,还包括比对的工具和对应的R包, featrueCounts的定量效果和HTSeq差不多,但是featureCounts的优点非常快!

4、cuffmerge整合所有组装好的转录本

当我们使用Cufflinks处理多个数据之后,我们需要将其转录本数据整合为一个全面的转录本集合,Cuffmerge是一个将Cufflinks生成的gtf文件融合为一个更加全面的转录本注释结果的工具。

cuffmerge -g annotation.gtf -s genome.fa -p 10 -o cuffmerge_output assembles.txt

-g  参考注释文件

-s  基因组文件

-o 输出文件夹

-p 使用的线程数

assembles.txt 提供组装的转录本位置

cat > assembles.txt

~/cufflink/C1_R1_output/transcripts.gtf

~/cufflink/C1_R2_output/transcripts.gtf

~/cufflink/C1_R3_output/transcripts.gtf

SAM文件

9778818威尼斯官网 6

20171218-215939.png

SAM代表Sequence Alignment/Map格式,是一种制表符分隔的文本格式,包含一个可选的头部分(header section,有人称之为“注释部分”),和一个比对部分(alignment section)。如果包含头部分,那么头部分必须置于比对部分之前。头部分的行以@符号开头,而比对部分的行不以@符号开头。比对部分的每一行包含11个必选的字段,用于说明重要的比对信息,如比对位置(mapping position)等;另有可变数量的可选字段,用于存储其他信息(flexible)或比对软件特异的信息。

(2)在Linux中通过验证MD5值检验文件是否完整

参考链接:https://www.cnblogs.com/zhuxiaohou110908/p/5786893.html

9778818威尼斯官网 7

featureCounts下载

二进制版本的软件解压即可使用

#下载subread
wget https://sourceforge.net/projects/subread/files/subread-1.5.3/subread-1.5.3-Linux-x86_64.tar.gz
tar zxvf subread-1.5.3-Linux-x86_64.tar.gz

5、cuffdiff计算差异表达的基因

计算的简单原理在于测序深度和外显子长度一定时,Read的数量与对应的转录本数量成正比。通过对Reads进行计数计算转录本的表达量。同时cuffdiff可以计算不同条件下转录本表达水平的显著性差异。

cuffdiff -o cuffdiff_output -p 10 -L Control_1,Control_2 -u cuffmerge_output/merged.gtf  ./tophat_output_C1_R1/accepted_hits.bam,./tophat_output_C1_R2/accepted_hits.bam ./tophat_output_C2_R1/accepted_hits.bam,./tophat_output_C2_R2/accepted_hits.bam

-p 使用的线程数

-u 同cuffmerge使用的-u参数

-L 为每个样品标上名称

接下来Cuffmerge产生的gtf文件,Cuffdiff需要它提供的注释进行初始转录产物和可变剪切等定量分析。最后是TopHat产生的bam文件,如果一个样品中有多个实验重复,那么我们需要提供哦呢bam文件列表,文件名之间以逗号隔开。

运行之后,cuffdiff输出的文件在diff_out目录之下。其中包括一些按类别统计的表达水平结果,如果有相同的转录起始位点,或具有相同的编码区的转录本的表达水平,我么你可以利用他们进行下一步的分析

header部分

@开头,后面紧接着一个或两个字母的头记录类型编码。每行都是以制表符分隔的,除了@CO行,每个数据字段都服从‘TAG:VALUE’这一格式,其中TAG是一个两字符的字符串,定义了VALUE的格式和内容。所以头部分的行都符合这一正则表达式: /^@[A-Z][A-Z](t[A-Za-z][A-Za-z0-9]:[ -~] ) $/或者 /^@COt.*/

9778818威尼斯官网 8

header.png

3.2 将数据链接到野生型、分离群体混池文件夹中

$ cd /myhome/MutMap_test/1.qualify_read/anyname(野生型亲本)

$ ln –s /fastq/ZZZ/ZZZ_L003_R1_001.fastq.gz p_3_1_sequence.txt.gz

$ ln –s /fastq/ZZZ/ZZZ_L003_R2_001.fastq.gz p_3_2_sequence.txt.gz

$ ln –s /fastq/ZZZ/ZZZ_L004_R1_001.fastq.gz p_4_1_sequence.txt.gz

$ ln –s /fastq/ZZZ/ZZZ_L004_R2_001.fastq.gz p_4_2_sequence.txt.gz

$ cd /myhome/MutMap_test/1.qualify_read/mybulk(突变体混池)

$ ln –s /fastq/XXX/XXX_L005_R1_001.fastq.gz p_5_1_sequence.txt.gz

$ ln –s /fastq/XXX/XXX_L005_R2_001.fastq.gz p_5_2_sequence.txt.gz

$ ln –s /fastq/XXX/XXX_L006_R1_001.fastq.gz p_6_1_sequence.txt.gz

$ ln –s /fastq/XXX/XXX_L006_R2_001.fastq.gz p_6_2_sequence.txt.gz

HiCUP由6个Perl脚本组成,分别如下:

常用参数说明:

-p 针对paired-end数据
-T 多线程数
-t 默认将exon作为一个feature
-g 默认是gene_id,从注释文件中提取Meta-features信息用于read count
-a 基因组注释文件,默认是gtf
-o 输出文件

(PS:存储不够了,后续选择两组数据做流程分析。)

echo "hisat2 started at $(date)"
#make index 
cd /pnas/fangxd_group/renyx/macaque/00ref
hisat2-build -p 8 Macaca_mulatta.Mmul_8.0.1.dna.toplevel.fa maca_index

#alignment
for i in {230..231}
do
hisat2 -t -x 00ref/maca_index -1 02clean_data/SRR4042$i_paired_clean_R1.fastq.gz  
-2 02clean_data/SRR4042$i_paired_clean_R2.fastq.gz 
 -S 03align_out/SRR4042$i.sam 
done

#covert to bam, and sort bam
for i in {230..231}
do
samtools view -Sb SRR4042$i.sam > SRR4042$i.bam
samtools sort SRR4042$i.bam -o SRR4042$i_sorted.bam
samtools index SRR4042$i_sorted.bam 
rm *sam
done 
echo "hisat2 finished at $(date)"

#featureCounts定量
cd /pnas/fangxd_group/renyx/macaque/
featureCounts=/pnas/fangxd_group/renyx/software/subread-1.5.3-Linux-x86_64/bin/featureCounts
$featureCounts -p -T 6 -t exon -g gene_id -a 00ref/Macaca_mulatta.Mmul_8.0.1.91.gtf 
-o 05featurecount_out/counts.txt  03align_out/hisat2_mapping/*sorted.bam

6、差异基因的定义

awk '( $10 <-1 && $12 < 0.05 ){print $0}' gene_exp.diff > up_regulate.txt(上调)

awk '( $10 > 1 && $12 < 0.05 ){print $0}' gene_exp.diff > down_regulate.txt(下调)

比对部分

每行表示一个片段(segment)的线性比对(linear alignment)。每行有11个必选字段,信息缺失字段用0或者“*”表示。

9778818威尼斯官网 9

比对字段.png

4.放置参考基因组fasta文件

#将参考基因组fasta格式文件放在下面的文件夹中

/myhome/MutMap_test/downloaded_fasta/IRGSP-1.0_genome.fasta/

#为fasta格式参考基因组创建.fai文件

samtools faidx IRGSP-1.0_genome.fasta

HiCUP Digester:确定reference上的酶切位点

Hisat2的比对结果:

SRR4042230和SRR4042231的比对率分别为88.02%和88.81%,总比对时间将近5小时。

9778818威尼斯官网 10

Hisat2 mapping

FLAG字段的理解

SAM 文件是比对结果的展示, 格式详细介绍见官方文档, 其中第二列为按位标记的 flag, 这是一种常用且高效的保存多个布尔特征值的方法, 可通过按位取 '与' 和取 '或' 快速提取对应的值

举个简单的例子: 在 SAM 格式中, 当 flag 为 1, 也即对应的二进制为 01 时, 表示该 read 有多个测序数据 (template having multiple segments in sequencing), 一般理解为有双端测序数据 (另一条没被过滤掉), 而 flag 为 2, 也即二进制 10 时, 表示这条 read 的多个片断都有比对结果 (each segment properly aligned according to the aligner), 通常理解为双端 reads 都比对上了, 那么就可以推断出 flag 为 3 时, 也即二进制的 11, 表示该 read 有另一端的 read 并且比对成功, 可以看到, 其实就是 01 加 10, 反之, 当我们看到二进制的 11, 就可以想到这表示这两位的状态均为 1 (通常即为真), 再查询定义就知道了具体的含义

这样就通过一个位来反映一个 True/ False 的状态值, 同时方便了人脑和电脑, 在 SAM flag 中, 目前共有 12 种不同状态, 列举如下:
(来源:http://not.farbox.com/post/sam_flag)

9778818威尼斯官网 11

20171219-124131.png

SAM Format提供了工具可快速查询各个value的含义。
将其按需求自由组合, 再结合使用 samtools 软件 view 命令的 -f (满足) -F (不满足) 和 -c (计数) 参数, 就可以实现各种过滤和统计了, 比如 samtools view -f 4 就会提取出未比对上的 reads, samtools view -c -F 4 就能统计出比对上的 reads 数, 也即比对率的分子部分

5. 配置common.fnc文件

HiCUP:为主程序,依次执行以下步骤

featureCounts结果产生两个文件:

hisat_counts.txt.summary包含一些基本的统计信息。

9778818威尼斯官网 12

hisat_counts.txt.summary

hisat_counts.txt
包含8列信息,分别是Geneid;Chr 染色体号,这里会显示多个数目;Start; End; Strand; Length; 第7列和第8列显示的是Counts数

9778818威尼斯官网 13

hisat_counts.txt


CIGAR字段

9778818威尼斯官网 14

20171219-125551.png

5.1编辑config.tx文件

在第12、15、31、52、115行做出如下修改:

#在第12行修改突变体混池名称

Key1_Bulked_sample_name="mybulk"

#在第15行修野生型亲本名称

Key1_My_cultivar_sample="anyname"

#第31行指定参考基因组的路径

Key2_Path_public_reference_FASTA="./downloaded_fasta/IRGSP1.0_genome.fasta/IRGSP-1.0_genome.fasta"

#在第52行配置BWA的线程数(实验室的服务器最大32,设置为20)

Key2_BWA_CPU=20

#第115行设置突变体混池的个体数

Key3_Individuals=20

其他设置一般不需要修改,具体设置参考MutMap Protocol

HiCUP Truncater:在reads上寻找酶切位点,并将reads切开

subread 比对 定量

#make subread index
echo "subread  started at $(date)"
buildindex=/pnas/fangxd_group/renyx/software/subread-1.5.3-Linux-x86_64/bin/subread-buildindex
cd /pnas/fangxd_group/renyx/macaque/00ref
$buildindex -o maca Macaca_mulatta.Mmul_8.0.1.dna.toplevel.fa

#alignment
subjunc=/pnas/fangxd_group/renyx/software/subread-1.5.3-Linux-x86_64/bin/subjunc
subjunc_maca_index=/pnas/fangxd_group/renyx/macaque/00ref/maca
cd /pnas/fangxd_group/renyx/macaque/
for i in {230..231}
do
$subjunc -T 5 -i $subjunc_maca_index -r 02clean_data/SRR4042$i_paired_clean_R1.fastq.gz -R 02clean_data/SRR4042$i_paired_clean_R2.fastq.gz -o 03align_out/subread_mapping/SRR4042$
done

#counts expre
featureCounts=/pnas/fangxd_group/renyx/software/subread-1.5.3-Linux-x86_64/bin/featureCounts
$featureCounts -p -T 6 -t exon -g gene_id -a 00ref/Macaca_mulatta.Mmul_8.0.1.91.gtf 
-o 05featurecount_out/subread_counts.txt  03align_out/subread_mapping/*_subjunc.bam
echo "subread finished at $(date)"

9778818威尼斯官网 15

image.png

9778818威尼斯官网 16

image.png

可以看出subjunc比对不到2小时,比hisat2快3个多小时。


第10及11列就是我们FASTQ文件中的序列及phred 33对应的ASCII字符

5.2 加载配置

#切换进入MutMap top目录下

$ cd /myhome/MutMap_test

#运行Bat_make_common.fnc.sh

$ ./Bat_make_common.fnc.sh

# 上面的命令会做出如下的操作

mv 1.qualify_read/mybulk 1.qualify_read/XXX

mv 1.qualify_read/anyname 1.qualify_read/ZZZ

mkdir -p 1.qualify_read/XXX/q30p90/sep_pair

mkdir -p 1.qualify_read/ZZZ/q30p90/sep_pair

HiCUP Mapper:将reads比对到参考基因组上,如果输入的是PEreads,则R1和R2分开单独比对到reference上。比对内部调用bowtie或bowtie2比对。这一步会利用到事先建好的bowtie2 index

提取counts

根据第1列是Geneid,第7,8列是counts数,用awk提取出geneID和counts。

awk -F 't' '{print $1,$7,$8}' OFS='t' hisat_counts.txt >hisat_matrix.out
awk -F 't' '{print $1,$7,$8}' OFS='t' subread_counts.txt >subread_matrix.out

6.分析数据

HiCUP Filter:结合HiCUP Digester生成的酶切位点文件,过滤掉常见的Hi-C artefacts,例如Dangling Ends等

参考资料:

一个植物转录组项目的实战
史上最快的转录组流程-subread
转录组定量-FEATURECOUNTS

6.1 对测序数据进行质控和过滤

HiCUP Deduplicator:移除(仅保留一处最佳比对) PCR重复

(1)对野生型亲本测序数据进行质控和过滤

# cd进入1.qualify_9778818威尼斯官网,read文件夹下

$ cd /myhome/MutMap_test/1.qualify_read

#对野生型亲本测序数据进行质控和过滤

$ nohup ./Run_all_Bats.sh 9 &

#查看日志文件nohup.log,自动追加最新工作结果并打印到屏幕

tail -f nohup.log

  1. HiCUP的使用

(2)对突变体混池测序数据进行质控和过滤

#对突变体混池测序数据进行质控和过滤

$ nohup ./Run_all_Bats.sh 0 &

#查看日志文件nohup.log,自动追加最新工作结果并打印到屏幕

tail -f nohup.log

3.1 先采用bowtie2-build对reference建立索引

6.2 构建参考基因组

/path/to/bowtie2-build reference.fasta reference

(1)Run “Run_all_Bats.sh”

将野生型亲本测序reads比对到参考基因组上,并替换SNP,构建自己的参考基因组:

#cd进入2.make_consensus文件夹下

$ cd /myhome/MutMap_test/2.make_consensus/

#运行程序

$nohup ./Run_all_Bats.sh &

3.2 采用hicup目录下的hicup_digester在reference上寻找酶切位点,生成酶切信息文件

(2)Run “Run_all_Bats.sh”

将野生型亲本的reads重新比对到新构建的参考基因组上,发现错配造成的SNP,用于下一步的排除和过滤:

#cd进入下面的文件夹

$ cd /myhome/MutMap_test/2.make_consensus/90.align_to_this_fasta/

#运行如下程序

$ nohup ./Run_all_Bats.sh &

/path/to/hicup_v0.6.1/hicup_digester

6.3突变体混池数据分析及作图

#切换到3.analysis文件夹下

$ cd /myhome/MutMap_test/3.analysis/

#运行如下命令,进行数据分析和作图

$ nohup ./Run_all_Bats.sh &

--re1 ^GATC,MboI

7.数据存储位置

--genome reference

(1)构建的野生型参考基因组

数据存放在2.make_consensus/50.make_consensus中,格式为fasta

--outdir /path/to/output/dir

(2)突变体混池和野生型亲本比对鉴定出的snp信息

存储在3.analysis/70.awk_custom中

reference.fasta

(3)SNP-index图

3.analysis/90.slidingwindow/pngs

3.3 再采用hicup主程序进行分析

(4)SNP低密度区域

每个窗口中的SNP数目小于10:

3.analysis/90.slidingwindow/pngs/Z1MZ2K/mut_index_X/mask10

有两种方式运行hicup,其一是将所有参数写到config文件中,可以先运行

8.使用samtools查看突变位点

使用方法:

samtools tview align.bam ref.fasta

samtools tview命令说明:

按?查看帮助,q或Esc退出

空格键向右移动一屏,Backspace键向左移动一屏

.意为和正链匹配,,意为和反向链匹配

大写字母意为,正链替换为该碱基;小写字母意为反向炼替换为该碱基

按g可跳转至指定的碱基,如:8:1256321

也可以使用IGV查看突变位点。

hicup --example

9. 突变位点的注释

Pepline的输出结果中并没有直接给出突变位点的注释信息,不利于突变位点的筛选,得到所有的SNP后,我们后续可以通过snpEFF、ANNOVAR、AnnTools等软件,对突变进行注释。

生成样例config文件,修改其中的参数,并运行

10. 结果展示

此处展示一张我用MutMap pipeline分析出的结果,找到了突变基因的连锁区间和候选突变位点。

9778818威尼斯官网 17

或者直接用命令行运行,如下:

/path/to/hicup_v0.6.1/hicup

--bowtie2 /path/to/bowtie2

--digest Digest_reference_MboI_None_18-36-52_07-08-2018.txt

--format Sanger

--index /path/to/reference/index

--outdir /path/to/output/dir

--threads 40

/path/to/reads*.fastq.gz

Ø xx_R1_2.hicup.sam

Ø xx_ R1_2.HiCUP_summary_report.html

前者是最终的sam文件,后者是全程汇总报告

HiCUP运行后hicup_truncater先生成xx.R1.trunc.fastq和xx.R2.trunc.fastq文件,它是按照酶切位点对Reads进行切除,因此得到的reads长短不一。完成后会统计Truncated reads数,Not Truncated reads数,以两者的占比,总的reads数,平均reads长度,并记录在hicup_truncater_summary_xx.txt文件中。两个fastq对应的svg图像中仅绘制了前两项的条形图。

hicup_mapper调用bowtie2-align-s进行比对处理,生成xx_ R1_2.pair.sam文件,同时也统计了太短而无法比对的reads数,唯一比对的reads数,多处比对的reads数,比对不上的reads数,成对的reads数,以及以上各项的占比,记录在hicup_mapper_summary_xx.txt文件中,顺带还画了R1和R2的比对结果条形图(xx_R*_trunc.fastq.mapper_barchart.svg)

hicup_filter对上步生成的sam文件进行过滤,分别将每种invalid ditag类型(包括same internal,same dangling ends,same circularised,re_ligation及其它)过滤掉的比对结果写入到hicup_filter_ditag_rejects_xx/目录下

过滤后可用于下步分析的sam文件为xx_R1_2.filt.sam

当然,对结果ditag结果也进行了统计,包括valid pairs,invalid pairs,same circularised,same fragment dangling ends,same fragment internal, re-ligation, contiguous sequence, wrong size以及total pairs等类型,记录在hicup_filter_summary_xx.txt文件中。同时也画了Di-tag长度分布svg图和各种类型的svg piechart图。

9778818威尼斯官网 18

hicup_deduplicator去重,并且计算了序列内的unique reads(<10kb和>10kb分别统计)以及序列间的unique reads,并画了相应的svg图

最终的结果为xx_R1_2.hicup.sam

注意,这个最终的sam文件中仅存放De-duplication后Unique Di-Tags的Read Pairs,即如下图中的read pairs(强调:是read pairs!)

9778818威尼斯官网 19

最后不仅给出了所有步骤reads数据变化的汇总(HiCUP_summary_report_xx.txt),同时还给给出了非常美观的html报告(xx_ R1_2.HiCUP_summary_report.html)!

另外,需要强调一点,hicup的结果如果想接到Lachesis做组装。在hicup生成的sam转bam的时候一定要用sort -n,按名字排序!如果用默认的参考序列位置排序,则Lacheis中会生成Cluster信息,但是没有Order信息!

PE reads总数据量45G,参考基因组440M。运行指定40个核,估算出的CPU最大消耗35个整,最大消耗内存6.34G,运行时间约7.5h,最大消耗是hicup_mapper比对。

9778818威尼斯官网 20

本文由9778818威尼斯官网发布于9778818威尼斯官网,转载请注明出处:HiCUP的应用以至结果解读,Seq数据深入分析

您可能还会对下面的文章感兴趣: