澳门新浦京娱乐场网站-www.146.net-新浦京娱乐场官网
做最好的网站

澳门新浦京娱乐场网站:转录组笔记,转录组入

任务列表

参考http://www.jianshu.com/p/681e02e7f9af
http://www.jianshu.com/p/93f96e7538da
任务:

  1. 比对软件很多,首先大家去收集一下,因为我们是带大家入门,请统一用hisat2,并且搞懂它的用法。
  2. 直接去hisat2的主页下载index文件即可,然后把fastq格式的reads比对上去得到sam文件。
  3. 接着用samtools把它转为bam文件,并且排序(注意N和P两种排序区别)索引好,载入IGV,再截图几个基因看看!
  4. 顺便对bam文件进行简单QC,参考直播我的基因组系列。

HISAT2:比对到基因组

  • 比对软件
  • hisat2的用法
  • 下载index文件
  • 比对、排序、索引
  • 质量控制
  • 载入IGV,截图几个基因

1. 比对软件

  • HISAT2:http://ccb.jhu.edu/software/hisat2/index.shtml
    参考资料:http://blog.biochen.com/archives/337
  • STAR:https://codeload.github.com/alexdobin/STAR/zip/master
    参考资料:http://www.bio-info-trainee.com/727.html
  • TopHat:http://ccb.jhu.edu/software/tophat/index.shtml
    参考资料:http://blog.sina.com.cn/s/blog_8808cae20101amqp.html
  • RapMap:https://github.com/COMBINE-lab/RapMap
    参考文献:https://academic.oup.com/bioinformatics/article/32/12/i192/2288985/RapMap-a-rapid-sensitive-and-accurate-tool-for
  • CIDANE:http://ccb.jhu.edu/software/cidane/
    参考文献:https://genomebiology.biomedcentral.com/articles/10.1186/s13059-015-0865-0
  • CLASS2 :https://sourceforge.net/projects/splicebox/files/?source=navbar
    参考文献:https://academic.oup.com/nar/article/44/10/e98/2516329/CLASS2-accurate-and-efficient-splice-variant

TopHat首次被发表已经是7年前,STAR的比对速度是TopHat的50倍,HISAT更是STAR的1.2倍。HISAT2是TopHat2/Bowti2的继任者,使用改进的BWT算法,实现了更快的速度和更少的资源占用,作者推荐TopHat2/Bowti2和HISAT的用户转换到HISAT2。
官网:https://ccb.jhu.edu/software/hisat2/index.shtml(学习一个软件最好的方法就是结合现有中文资料,加上阅读官方说明书和HELP文档,一般刚开始学习的时候先使用默认参数,不要乱调参数)

hisat2的用法

2. HISAT2的使用

为什么要用index?官网有描述:为了用整个index代表整个基因组,HISAT2 用小的index覆盖了整个基因组,每个index覆盖了56 Kbp的范围,覆盖整个人类基因组需要55,000 indexes,这些index结合其他策略可以快速准确的比对序列。

#index 下载
nohup wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz &
tar -zxvf *.tar.gz #解压
# 删除压缩包
rm -rf *.tar.gz
  • hisat2 -h查看帮助文件:
Usage: 
  hisat [options]* -x <bt2-idx> {-1 <m1> -2 <m2> | -U <r> | --sra-acc <SRA accession number>} [-S <sam>]

  <bt2-idx>  Index filename prefix (minus trailing .X.ht2).
  <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).
  <SRA accession number>        Comma-separated list of SRA accession numbers, e.g. --sra-acc SRR353653,SRR353654.
  <sam>      File for SAM output (default: stdout)

  <m1>, <m2>, <r> can be comma-separated lists (no whitespace) and can be  specified many times.  E.g. '-U file1.fq,file2.fq -U file3.fq'.
  • 参数说明:

-x 指定index文件名
-1 <m1> 双端测序第一个文件
-2 <m2> 双端测序第二个文件
-U 单端测序文件
--sra-acc 登录号
-S 输出文件为sam格式
-q 输入文件为FASTQ .fq/.fastq格式

-比对到参考基因组,RNA-Seq数据是从 SRR3589957 ~ SRR3589962,6个样本的PE数据,也就是有6次循环,可以写脚本,也可以直接在命令行里运行

for i in `seq 56 62`
do
    hisat2 -t -x /opt/NfsDir/UserDir/qin/qin/Data/annotation/hg19/genome -1 /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/SRR35899${i}.sra_1.fastq.gz -2 /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/SRR35899${i}.sra_2.fastq.gz -S /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/SRR35899${i}.sam &
done

运行时间如下:

澳门新浦京娱乐场网站 1

Paste_Image.png

  • 很耗CPU啊!用的服务器!

澳门新浦京娱乐场网站 2

Paste_Image.png

  • 结果:

澳门新浦京娱乐场网站 3

Paste_Image.png

SAM(sequence Alignment/mapping)数据格式是目前高通量测序中存放比对数据的标准格式,当然他可以用于存放未比对的数据。目前处理SAM格式的工具主要是SAMTools。samtools功能众多,在本次作业中,我们主要学会将sam文件转换为bam文件,并对bam文件进行sorted(其中有两种排序方式N和P),最后建立索引。

Program: samtools (Tools for alignments in the SAM format)
Version: 1.3.1-58-gcbee45e (using htslib 1.3.2-228-g0c32631)

Usage:   samtools <command> [options]

Commands:
  -- Indexing
     dict           create a sequence dictionary file
     faidx          index/extract FASTA
     index          index alignment

  -- Editing
     calmd          recalculate MD/NM tags and '=' bases
     fixmate        fix mate information
     reheader       replace BAM header
     rmdup          remove PCR duplicates #移除PCR产生的重复序列
     targetcut      cut fosmid regions (for fosmid pool only)
     addreplacerg   adds or replaces RG tags

  -- File operations
     collate        shuffle and group alignments by name
     cat            concatenate BAMs
     merge          merge sorted alignments
     mpileup        multi-way pileup
     sort           sort alignment file
     split          splits a file by read group
     quickcheck     quickly check if SAM/BAM/CRAM file appears intact
     fastq          converts a BAM to a FASTQ #格式转换
     fasta          converts a BAM to a FASTA

  -- Statistics
     bedcov         read depth per BED region #bed文件的测序深度
     depth          compute the depth 
     flagstat       simple stats
     idxstats       BAM index stats
     phase          phase heterozygotes
     stats          generate stats (former bamcheck)

  -- Viewing
     flags          explain BAM flags
     tview          text alignment viewer
     view           SAM<->BAM<->CRAM conversion
     depad          convert padded BAM to unpadded BAM

主要功能:
view: BAM-SAM/SAM-BAM 转换和提取部分比对
sort: 比对排序
merge: 聚合多个排序比对
index: 索引排序比对
faidx: 建立FASTA索引,提取部分序列
tview: 文本格式查看序列
pileup: 产生基于位置的结果和 consensus/indel calling

# 首先将比对后的sam文件转换成bam文件
# 利用的是samtools的view选项,参数-S 输入sam文件;参数-b 指定输出的文件为bam;最后重定向写入bam文件
$ for ((i=56;i<=62;i  ));do samtools view -S /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/SRR35899${i}.sam -b > /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/SRR35899${i}.bam;done
# 将所有的bam文件进行排序
$ nohup for (( i=58;i<=62;i   )); do samtools sort /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/SRR35899${i}.bam -o /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/SRR35899${i}.sorted.bam;done
# 将所有的排序文件建立索引,索引文件.bai后缀
$ for ((i=56;i<=62;i  ));do samtools index /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/SRR35899${i}.sorted.bam;done

合在一块跑:
for i in `seq 56 58`
do
    samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam
    samtools sort SRR35899${i}.bam -o SRR35899${i}_sorted.bam
    samtools index SRR35899${i}_sorted.bam
done

澳门新浦京娱乐场网站 4

Paste_Image.png

直接去hisat2的主页下载index文件即可,然后把fastq格式的reads比对上去得到sam文件。 接着用samtools把它转为bam文件,并且排序(注意N和P两种排序区别)索引好,载入IGV,再截图几个基因看看! 顺便对bam文件进行简单QC

本作业是比对到基因组,所以使用gapped or splices mapper,此流程已经更新。TopHat首次被发表已经是7年前,STAR的比对速度是TopHat的50倍,HISAT更是STAR的1.2倍。HISAT2是TopHat2/Bowti2的继任者,使用改进的BWT算法,实现了更快的速度和更少的资源占用,作者推荐TopHat2/Bowti2和HISAT的用户转换到HISAT2。

比对质控(QC)

常用工具有

  • Picard https://broadinstitute.github.io/picard/
  • RSeQC http://rseqc.sourceforge.net/
  • Qualimap http://qualimap.bioinfo.cipf.es/
java -jar /opt/NfsDir/BioDir/picard-tools-1.119/picard.jar CollectMultipleMetrics 
      I=/opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/SRR3589956.sorted.bam 
      O=/opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/multiple_metrics 
      R=/opt/NfsDir/PublicDir/reference/ucsc.hg19.fasta 

统计bam文件:

/opt/NfsDir/BioDir/RSeQC/RSeQC-2.6.4/scripts/bam_stat.py -i /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/SRR3589956.sorted.bam

提示出错:

澳门新浦京娱乐场网站 5

Paste_Image.png

检查发现是路径和名称写错了,修改后就可以了
对bam文件进行统计分析

澳门新浦京娱乐场网站 6

Paste_Image.png

#下载hg19_RefSeq.bed文件
https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_RefSeq.bed.gz/download
#查看基因组覆盖率
/opt/NfsDir/BioDir/RSeQC/RSeQC-2.6.4/scripts/read_distribution.py -i /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/SRR3589956.sorted.bam -r /opt/NfsDir/UserDir/qin/qin/Data/annotation/hg19/hg19_RefSeq.bed

澳门新浦京娱乐场网站 7

Paste_Image.png

  1. 下载index文件(小鼠)
    axel ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/mm10.tar.gz
    tar -xvzf mm10.tar.gz

    下载注释文件
    axel ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_mouse/release_澳门新浦京娱乐场网站,澳门新浦京娱乐场网站:转录组笔记,转录组入门。M10/gencode.vM10.annotation.gtf.gz

  2. 序列比对
    a. 配置 1核8G内存
    hisat2 -t -x ref/mm10/genome -1 rawdata/SRR3589960.sra_1.fastq.gz -2 rawdata/SRR3589960.sra_2.fastq.gz -S align/SRR3589960.sam

官网:

![](https://upload-images.jianshu.io/upload_images/10354448-7e1d4e5f0a826b5d.jpg)

111111111111.jpg



运行时间38分钟,源文件1.6G

下载index文件

更改配置4核32G
for ((i=60;i<=62;i=i ));do hisat2 -t -x ref/mm10/genome -1 rawdata/SRR35899${i}.sra_1.fastq.gz -2 rawdata/SRR35899${i}.sra_2.fastq.gz -S align/SRR35899${i}.sam;done

cd ~/reference
mkdir -p index/hisat && cd index/hisat
wget -c ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz
wget -c ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/mm10.tar.gz
tar zxvf hg19.tar.gz
tar xvzf mm10.tar.gz

注释: -t 记录时间 -x hg19(index)文件路径 -1 -2 测序的两个fastq文件 -S 比对结果输出路径 -U 单端测序文件
reference=~/wikiwei/human/ref/index/hg19/genome
hisat2 -t -x $reference -U SRR957679.fastq -S siSUZ12_1.sam 2>siSUZ12_1.log

-c:断点续传

参考文章

比对、排序、索引

  1. 转录组入门(mac版本)
  2. 澳门新浦京娱乐场网站:转录组笔记,转录组入门。https://www.cnblogs.com/freescience/p/7342895.html

把fastq格式的reads比对上去得到sam文件,接着用samtools把它转为bam文件,并且排序(注意N和P两种排序区别)索引好(可以使用管道实现,省去中间SAM保存的过程,直接输出BAM文件)

编写bash脚本:map.sh

#! usr/bin/bash
set -u
set -e
set -o pipefail
hg19_ref=/mnt/hgfs/2017/reference/index/hisat/hg19/genome
mm10_ref=/mnt/hgfs/2017/reference/index/hisat/mm10/genome
data_path=/mnt/hgfs/2017/rna_seq/data
NUM_THREADS=25
ls --color=never Homo*1.fastq.gz | while read id;do(~/biosoft/hisat2-2.1.0/hisat2 -t -p $NUM_THREADS -x $hg19_ref -1 $data_path/${id%_*}_1.fastq.gz -2 $data_path/${id%_*}_2.fastq.gz 2 > ${id%_*}_map.log | samtools view -Sb  - > ${id%_*}.bam);done
ls --color=never Mus*1.fastq.gz | while read id;do(~/biosoft/hisat2-2.1.0/hisat2 -t -p $NUM_THREADS -x $mm10_ref -1 $data_path/${id%_*}_1.fastq.gz -2 $data_path/${id%_*}_2.fastq.gz 2 > ${id%_*}_map.log | samtools view -Sb  - > ${id%_*}.bam);done
ls --color=never *.bam | while read id;do(samtools sort --threads $NUM_THREADS $id -o ${id%.*}_sorted.bam);done
ls --color=never *_sorted.bam | while read id;do(samtools index $id);done

运行脚本: 

bash map.sh

质量控制

对bam文件进行简单QC

Reads比对后的质量控制(评估比对质量的指标)

比对上的reads占总reads的百分比;

Reads比对到外显子和参考链上的覆盖度是否一致;

比对到基因组序列,多重比对reads;

相关质控软件除了Picard,RSeQC,Qualimap还有一大堆

本文由澳门新浦京娱乐场网站发布于澳门新浦京娱乐场网站,转载请注明出处:澳门新浦京娱乐场网站:转录组笔记,转录组入