想用smrt analysis-2.3的过来看看

最近找smrt analysis 2.3的程序,真是太辛苦了。我不想用SMRT的图形界面,pacbio又把github上的相关项目obsolete了,tofu项目也没了,pbtranscript装起来各种错误,实在要放弃了,还好看到有人做了2.3的docker镜像。用最近的SMRT6.0中的isoseq3的同学可以忽略本文,isoseq3请移步https://github.com/PacificBiosciences/IsoSeq3/blob/master/README_v3.1.md

docker pull bowhan/smrtanalysis-2.3.0.140936
docker run -it -v /home/bowhan/myData:/data bowhan/smrtanalysis-2.3.0.140936 bash
# -v 挂载目录

#进入smrtshell,有smrtanalysis的环境
/opt/smrtanalysis/current/smrtcmds/bin/smrtshell
# 以iso-seq数据处理为例
# 1. CCS 根据fofn中的h5文件,生成CCS
ConsensusTools.sh CircularConsensus  \
    --minFullPasses 0  --minPredictedAccuracy 75 \
    --parameters /opt/smrtanalysis/install/smrtanalysis_2.3.0.140936/analysis/etc/algorithm_parameters/2014-09/ \
    --numThreads 8 --fofn input.fofn \
    -o /data/output

# 2. classify
pbtranscript.py classify \
    --cpus 8 --min_seq_len 300 \
    --flnc ${prefix}.isoseq_flnc.fasta \
    --nfl ${prefix}.isoseq_nfl.fasta \
    -d ${prefix}.out \
    ${prefix}.ccs.fa \
    ${prefix}.isoseq_draft.fasta 

# 3. cluster
pbtranscript.py cluster \
    isoseq_flnc.fasta final.consensus.fa \
    --nfl_fa isoseq_nfl.fasta -d clusterOut \
    --ccs_fofn reads_of_insert.fofn --bas_fofn input.fofn \
    --cDNA_size under1k --quiver --blasr_nproc 8 \
    --quiver_nproc 8 

本文转自https://hub.docker.com/r/bowhan/smrtanalysis-2.3.0.140936/
其他参考:https://github.com/ben-lerch/IsoSeq-3.0

一步到位下载hg19基因组文件

hg19对应GRCh37,UCSC提供hg19的参考基因组下载。UCSC的下载地址在ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/

需要经过下载每个染色体,然后解压合并成一个整个的基因组文件
ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/

其实这样有点浪费时间,还要考虑合并的时候染色体的顺序是否按照1,2,3而不是1,10,11排下来的。目前我知道的最简单的办法的,从GATK bundle中下载。比如hg19整个基因组的文件。下面是一步到位的命令,包括了fasta,fai,dict文件。

wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/ucsc.hg19*

GATK bundle还提供一下其他文件,可以看看( ftp://ftp.broadinstitute.org/bundle ),比如dict文件,hg38文件等。当然构建参考基因组不一定非要合并染色体,个人习惯。

补充:UCSC也提供一个genome.fa的下载,ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

Generate Accurate Consensus Sequences from a Single SMRTbell

Circular-Consensus-Sequence

bax2bam

bax2bam -o mynewbam mydata.1.bax.h5 mydata.2.bax.h5 mydata.3.bax.h5

circus consensus sequence

Circular-Consensus-Sequence

ccs –minLength=100 myData.subreads.bam myResult.bam

进行CCS由pacbio测序系统决定的,插入序列测多遍之后,可以用来校正随机错误。

参考

https://github.com/PacificBiosciences/unanimity/blob/develop/doc/PBCCS.md

SpliceMap官网示例教程

SpliceMap是一个从头开始发现和比对splice junction的工具。它提供高敏感度并且支持任意长度RNA-seq 序列片段read长度. SpliceMap将RNA-seq reads比对到参考基因组上用于发现splicing junctions. 它至少拥有与当前技术条件下其它分析工具同等的灵敏度和特异性。
官网 http://web.stanford.edu/group/wonglab/SpliceMap/manual.html
下载 http://web.stanford.edu/group/wonglab/SpliceMap/download.html
示例教程 https://web.stanford.edu/group/wonglab/SpliceMap/tutorial.html
以官网SpliceMap 3.3.5.2 example (Linux-x86 64bit)为例,介绍如何用来自21号染色体的100k 条100bp的RNA reads寻找junction。

测试是否正常./bin/runSpliceMap,如果报错的话,需要到src文件夹运行 ./install.sh ../bin来安装SpliceMap,运行./install-bowtie.sh ../bin来安装

示例文件夹下面的结构

ls SpliceMap3352_example_linux-64
all.gene.refFlat.txt  bin  data  genome  INSTALL  LICENSE  output  run.cfg  src  temp

run.cfg

最重要的文件,为文本文件,包含测序reads路径,基因组文件路径,配置设置等详见
https://web.stanford.edu/group/wonglab/SpliceMap/cfg.html,需要对每个数据集进行配置。在cfg文件中,注释以#开头,赋值用tag = value表示,比如genome_dir = ../../genome/hg18,多值用如下表示

> tag
value_1
value_2
...
<

genome directory

该路径下包括物种的染色体,bowtie的索引,测试示例中仅放了21号染色体和相关bowtie的索引。

ls genome                                                      
chr21.1.ebwt  chr21.2.ebwt  chr21.3.ebwt  chr21.4.ebwt  chr21.fa  chr21.rev.1.ebwt  chr21.rev.2.ebwt

data directory

样本的测序reads,真实使用中可以是其他位置。

temp directory

在SpliceMap运行过程中的临时文件,最初的序列比对结果存储在这里,所以这个文件夹可能会很大。

output directory

SpliceMap的输出结果文件夹。

all.gene.refFlat.txt

包含来自ensembl,refseq,knowGene数据所有基因的注释,可以用来寻找新的junction。

bin

SpliceMap的运行程序,不需要安装,可以拷贝到其他地方。

src

包含SpliceMap/Bowtie源程序的文件夹。

运行SpliceMap

在example文件夹中,运行大概两到三分钟可以结束

./bin/runSpliceMap run.cfg

输出结果文件夹output

junction_color.bed

包含发现的junctions,novel junction在bed文件中用红色表示,可以用UCSC基因组浏览器或cisGenome浏览器看。

track   name=junctions  description="SpliceMap junctions" itemRgb="On"
chr21   14669961        14672440        (2)[49_2](2/0)  50      -       14669961        14672440        96,96,96        2 50,50    0,2429
chr21   16125769        16127587        (1)[1_1](1/0)   50      +       16125769        16127587        192,192,192     2 50,50    0,1768

junction_color.new.bed

不在all.gene.refFlat.txt文件中的新的junction。

track   name=junctions  description="SpliceMap junctions" itemRgb="On"
chr21   17186728        17208745        (2)[6_2](2/0)   50      +       17186728        17208745        255,0,0 2       50,50      0,21967
chr21   23658336        23666569        (1)[1_1](1/0)   50      -       23658336        23666569        255,192,192     2 50,50    0,8183

junction_nUM_color.bed

junction_color.bed文件中的至少被一条唯一比对的read支持的junction。

track   name=junctions  description="SpliceMap junctions" itemRgb="On"
chr21   14669961        14672440        (2)[49_2](2/0)  50      -       14669961        14672440        96,96,96        2 50,50    0,2429
chr21   16125769        16127587        (1)[1_1](1/0)   50      +       16125769        16127587        192,192,192     2 50,50    0,1768

junction_nUM_color.new.bed

在junction_nUM_color.bed中且不在all.gene.refFlat.txt中的junction。

track   name=junctions  description="SpliceMap junctions" itemRgb="On"
chr21   17186728        17208745        (2)[6_2](2/0)   50      +       17186728        17208745        255,0,0 2       50,50      0,21967
chr21   23658336        23666569        (1)[1_1](1/0)   50      -       23658336        23666569        255,192,192     2 50,50    0,8183

coverage_up.wig

染色体位置中唯一比对的reads的覆盖度

track   type=bedGraph   name="SpliceMap unique coverage"
chr21   9944337 9944421 1
chr21   14070536        14070619        1

coverage_down.wig

染色体位置中多重比对的reads的覆盖度

track   type=bedGraph   name="SpliceMap non-unique coverage"
chr21   16276385        16276484        -0.5
chr21   16277073        16277137        -0.5

coverage_all.wig

所以比对上的reads在染色体上的覆盖度

track   type=bedGraph   name="SpliceMap all coverage"
chr21   9944337 9944421 1
chr21   14070536        14070619        1

good_hits.sam

SAM格式存储的比对好的reads

R[87976]        73      chr21   9944337 255     85M     *       0       0       CATGTATCAACACCAAGTGCCATTTCTCTTAAGGATGCAAGGATGTTTTTTCCAAGATGGCAGATTACAGGCTTTTAGCATGCCT      IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII      NM:i:0  XC:i:15
R[28405]        177     chr21   14070536        255     100M    =       14070620        -15     GCAGAGGCTATTATTCATGGACTATCCAGTCTAACAGCTTGCCAGTTGAGAACGATATACATATGTCAGTTTTTGACAAGAATTTCAGCAGGAAAAACCC       IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII       MD:Z:40A43G15   NM:i:2  XC:i:0
R[28405]        113     chr21   14070620        255     100M    =       14070536        -15     GCAGCAGGAAAAACCCTTGATGCACAGTTTGAAAATGATGAACGAATTACACCCTTGGAATCAGCCCTGATGATTTGGGGTTCAATTGAAAAGGAACATG       IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII       MD:Z:42T42C14   NM:i:2  XC:i:0

debug_logs

输出log

用自己的数据

重点修改cfg文件中genome_dir (参考基因组路径),reads_list1和2(fastq文件),annotations (基因注释文件),temp_path(中间数据存放路径),out_path(结果文件夹),bowtie_base_dir (bowtie2的index文件,如果没有,会在temp_path生成)

参考

https://web.stanford.edu/group/wonglab/SpliceMap/tutorial.html
大部分翻译这个tutorial,具体数据是自己跑tutorial得到的。

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

Pacbio三代测序Primary Analysis Data文件夹

三代测序很多年了,刚工作的时候在超算中心做过三代的拼接,没好好研究过之后就再也没接触过,现在要做三代的项目,从头学习,Primary Analysis Data为初步数据分析文件夹,类似下面的文件夹结构

/path/to/secondary/storage/2420294/0011
├── Analysis_Results
│   ├── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.1.bax.h5
│   ├── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.1.log
│   ├── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.1.subreads.fasta
│   ├── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.1.subreads.fastq
│   ├── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.2.bax.h5
│   ├── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.2.log
│   ├── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.2.subreads.fasta
│   ├── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.2.subreads.fastq
│   ├── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.3.bax.h5
│   ├── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.3.log
│   ├── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.3.subreads.fasta
│   ├── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.3.subreads.fastq
│   ├── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.bas.h5
│   ├── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.sts.csv
│   └── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.sts.xml
├── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.1.xfer.xml
├── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.2.xfer.xml
├── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.3.xfer.xml
├── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.mcd.h5
└── m140415_143853_42175_c100635972550000001823121909121417_s1_p0.metadata.xml

主要文件有

bas.h5文件和bax.h5文件

bas.h5和相关的bax.h5文件是PacBio@RS II初级分析(primary analysis)的主要输出文件,这些文件由设备产生到本地存储位置,作为后续SMRT分析软件进行alignment、consensus和variant分析的输入文件。
PacBio@RS II之前,单个bas.h5文件包含了所有测序数据,随着PacBio@RS II升级,通量和read长度都在增加,现在包含一个bas.h5和3个bax.h5文件(1-3.bax.h5)。bax.h5文件包含测序的base call的信息,bas.h5是三个bax.h5的重要指针。
用h5dupm -n [movie name].bas.h5命令看一下文件

FILE_CONTENTS {
 group      /
 group      /MultiPart
 dataset    /MultiPart/HoleLookup
 dataset    /MultiPart/Parts
 }

可以看到有一个MultiPart的group组,下面有两个数据集,一个是HoleLookup,一个是Parts。
/MultiPart/Parts提供相关对应的的三个bax.h5的文件名
/MultiPart/HoleLookup对应ZMW零模波导孔编号

单个的bax.h5文件结构和原来一致 ,有两个组,分别是PulseData和ScanData。ScanData
包含了用于测试的仪器信息,一般不会用到。PulseData下面重点看/PulseData/BaseCalls组
用命令h5dump -d /PulseData/BaseCalls/Basecall [movie name].bax.h5 | head -10 查看数据集basecall如下序列

DATASET "/PulseData/BaseCalls/Basecall" {
   DATATYPE  H5T_STD_U8LE
   DATASPACE  SIMPLE { ( 305612594 ) / ( H5S_UNLIMITED ) }
   DATA {
   (0): 84, 84, 67, 84, 84, 67, 84, 84, 84, 84, 67, 71, 84, 71, 71, 84, 84,
   (17): 84, 84, 71, 84, 67, 84, 84, 84, 67, 84, 71, 84, 84, 84, 67, 84, 84,
   (34): 84, 67, 84, 84, 71, 84, 84, 67, 84, 71, 67, 71, 71, 84, 84, 84, 84,
   (51): 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84, 84,
   (68): 84, 84, 84, 84, 71, 84, 84, 71, 84, 84, 84, 84, 84, 84, 84, 84, 84,

其中A=> 65, C=> 67, G=>71, T=>84,与ASCII编码一致。

用h5dump -d /PulseData/BaseCalls/QualityValue bax.h5 | head -10 查看数据集qualityvalue如下

DATASET "/PulseData/BaseCalls/QualityValue" {
   DATATYPE  H5T_STD_U8LE
   DATASPACE  SIMPLE { ( 305612594 ) / ( H5S_UNLIMITED ) }
   DATA {
   (0): 10, 11, 21, 12, 11, 9, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 14, 0, 0, 0, 10, 0, 1,
   (23): 0, 0, 15, 1, 1, 1, 0, 1, 1, 1, 1, 0, 2, 10, 0, 0, 0, 0, 0, 0, 14, 20, 1,
   (46): 1, 30, 1, 0, 0, 1, 22, 0, 19, 0, 1, 0, 0, 21, 1, 1, 1, 0, 0, 17, 0, 0, 1,
   (69): 1, 1, 30, 2, 1, 0, 2, 1, 0, 10, 0, 1, 1, 0, 1, 1, 13, 0, 0, 0, 1, 0, 0,
   (92): 1, 0, 5, 1, 0, 0, 1, 0, 2, 1, 0, 0, 21, 0, 1, 0, 0, 20, 0, 0, 0, 16, 1,

表示测错的概率,用的是phred socre,和Illumina类似。当然,我看的这个文件可以看出测序质量并不是很好。
具体的组和数据集的说明,可以参考:https://s3.amazonaws.com/files.pacb.com/software/instrument/2.0.0/bas.h5+Reference+Guide.pdf

metadata.xml

包含测序仪ID,样本名,测序酶,试剂等元数据。
参考文件:https://s3.amazonaws.com/files.pacb.com/software/instrument/2.0.0/Metadata+Output+Guide.pdf

sts.xml

包含了单个movie的概括统计
参考文件:https://s3.amazonaws.com/files.pacb.com/software/instrument/1.3.1/Statistics+Output+Guide.pdf

用h5dump查看pacbio数据

查看h5格式的文件可以用hdf5下的h5dump工具查看,安装如下

wget -c https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.10/hdf5-1.10.5/src/hdf5-1.10.5.tar.gz
tar -xzvf hdf5-1.10.5.tar.gz
cd hdf5-1.10.5  
./configure &&make && make install

参考:

https://www.cnblogs.com/xudongliang/p/6908953.html
https://smrt-analysis.readthedocs.io/en/latest/Data-files-you-received-from-your-service-provider/

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

本条目发布于2019年3月20日。属于Default Category分类,被贴了 Pacbio、Sequencing 标签