# 01三代测序Iso-seq ## PacBio三代全长转录组测序 分析流程 ![分析流程](https://github.com/Magdoll/images_public/raw/master/github_isoseq3_wiki_figures/IsoSeq3_workflow_v3.png) 这张图片形象的展示了每一步测序数据的变化 ![测序说明](https://user-images.githubusercontent.com/39325949/66532093-b8a89100-eb40-11e9-8163-4a382a1afb53.png) ## 1.安装SMART软件 从网站[https://www.pacb.com/support/software-downloads/下载SMART-Link软件](https://www.pacb.com/support/software-downloads/%E4%B8%8B%E8%BD%BDSMART-Link%E8%BD%AF%E4%BB%B6) ```bash ./smrtlink_7.0.1.66975.run --rootdir /public/home/zpliu/software/smrtlink --smrttools-only ``` * \--rootdir 指定安装路径 * \--smrttools-only 只安装命令行工具，反正服务器里你图形界面也看不到 版本升级，下载升级版本的安装文件，然后只需要在安装的命令上面加个参数就行 ```bash ./新版本的安装程序 --rootdir /public/home/zpliu/software/smrtlink --smrttools-only --upgrade ``` 之后将安装的命令添加到环境变量中 ```bash export PATH=" /public/home/zpliu/software/smrtlink/smrtcmds/bin:$PATH" ``` ## 2.使用**CCS**对原始数据进行过滤 ```bash ccs --noPolish --minPasses 1 --minLength 300 --minSnr 4 --maxLength 10000 --maxDropFraction 0.8 --minPredictedAccuracy 0.8 --numThreads 20 --logFile ccs.log --reportFile repoet_css.txt input_subreads.bam output.bam ``` * **noPolish** 不会对数据进一步的过滤 * **-minPasses** 最低的通过值 * **minLength** 获取的draft consensus最短长度，用于下一步的分析 * **maxLength** 最长的长度 * **minSnr** 移除包含delete的SNP * **-minPredictedAccuracy** 最小精度0.8 * **--logFile** 记录日志文件 * **-reportFile** 报告处理的文件 * **maxDropFraction** Maximum fraction of subreads dropped by polishing ## 3.对转录本进行无参考基因组的归类 引物文件是固定的 ```bash $ cat primers.fasta >primer_5p AAGCAGTGGTATCAACGCAGAGTACATGGGG >primer_3p AAGCAGTGGTATCAACGCAGAGTAC ``` 3.1 去除引物 `lima output.bam primers.fasta demux.ccs.bam --isoseq --no-pbi -j 线程数 --min-length 300` 3.2 去除full length 的噪音remove polyA tails ```bash #v7版本的命令 isoseq3 refine 前缀primer_5p--primer_3p.bam primers.fasta movie.flnc.bam #v8版 isoseq3 refine movie.fl.P5--P3.bam primers.fasta movie.flnc.bam ``` 3.3 聚类 \` isoseq3 cluster movie.flnc.bam unpolished.bam · :warning:~~如果想要比较同源基因之间的差异的话，这一步可以不做~~ 3.4打磨 polich `isoseq3 polish -j 20 unpolished.bam input_subreads.bam polished.bam` ## 4.比对到参考基因组 4.1 软件安装 ```bash #下载最新版软件 ./configure --prefix=自定义安装路径 make make install ``` 4.2 建立参考基因组 ```bash gmap_build -D 存放索引的目录 -d G.arboreum.Chr.v1.0 基因组fasta文件 ``` 4.3 将全长转录本比对到参考基因组 ```bash gmap -D 索引所在目录 -d 索引文件前缀 -f samse -t 10 -n 2 polished.hq.fasta >gmap.sam 2>gmap.err ## 推荐使用这种方法,不然老是报错 cat fasta|gmap -D 索引所在目录 -d 索引文件前缀 -f samse -t 10 -n 2 >gmap.sam 2>gmap.err ``` * -f 输出文件为sam格式 * -t 指定线程数目 * -n 设置比对的类型，为0可以鉴定嵌合基因 4.4 根据SNP数据区分多倍体reads ```bash transcript/963 16 Chr01 1836618 0 62S165M1I453M74N371M119N348M172N * 0 0 GATCTCTACCATAAGCTTTTAGCAATGCCAAAATCAGTAAGATGAGCCTCTAAATCCTTGTCCAACAAAATATTTGATGACTTCACATCCCTGTGGATGATTCGTGGACTACAGTCATGATGTAAATACGCTAACCCTTGTGCAGCTCCCAGTGCAATCTTTAATCGAATGTTCCAACTGAGAACTTTCTTCTTGGTAGAAACATGGAGGAGATCCCAGAGACTGCCATTTTCCATATAGTCATAGAAGAGAAGGTTCCGAGACGGGGAGAGAGAATACCCTTGAAGGCTGACCAGATTTCGGTGCTTAATACTCCCAATTGTCTCGAGTTCTGTCTCGAATTCCTTCAAGCATTGTGGATAGTGAGAGTAGAGCCACTTGATGGCAACTGGCCT * MD:Z:97T72C38A23T42G109G36C1 5C47A89G20T100G16G141G77A11C10C5C10C17 NH:i:1 HI:i:1 NM:i:62 SM:i:40 XQ:i:24 X2:i:0 XO:Z:UU XS:A:- ``` 对于mapping状态可分为以下几类： * M：alignment match (can be a sequence match or mismatch) 表示read可mapping到第三列的序列上，则read的碱基序列与第三列的序列碱基相同，表示正常的mapping结果，M表示完全匹配，但是无论reads与序列的正确匹配或是错误匹配该位置都显示为M * I：insertion to the reference 表示read的碱基序列相对于第三列的RNAME序列，有碱基的插入 * D：deletion from the reference 表示read的碱基序列相对于第三列的RNAME序列，有碱基的删除 * N：skipped region from the reference 表示可变剪接位置 * P：padding (silent deletion from padded reference) * S：soft clipping (clipped sequences present in SEQ) * H：hard clipping (clipped sequences NOT present in SEQ) clipped均表示一条read的序列被分开，之所以被分开，是因为read的一部分序列能匹配到第三列的RNAME序列上，而被分开的那部分不能匹配到RNAME序列上。 * "="表示正确匹配到序列上 * "X"表示错误匹配到序列上 ## 5.TaMa将很相似的转录本合并**去冗余** ![合并转录本](https://github.com/GenomeRIK/tama/raw/master/images/Collapse1.png) 这个过程很复杂，图中就有两种可能的合并方式： * Transcription Start Site Collapse * Exon Cascade Collapse 具体可以查看这篇文献 5.1 下载和安装 ```bash #发现我没有pip2安装走一波 #下载文件 wget https://bootstrap.pypa.io/get-pip.py --no-check-certificate #执行安装 python get-pip.py --user git clone https://github.com/GenomeRIK/tama #也可以之间在windows下载好rz 进去 #程序基于python2运行 ，需要安装Biopython module模块 pip install --upgrade pip --user # 更新pip pip install Biopython --user ``` 运行tama\_collapse.py 脚本 ```bash /usr/bin/python tama_collapse.py -s ../gmap_sort.sam -f ../../Gr_genome/Graimondii_221_v2.0.fa -p tama -x capped ``` 每个参数的详细说明 ## 6.Cupcake去除冗余，这个步骤和5是一样的 **推荐这个流程** 6.1 安装Cupcake软件，这个流程适合依赖于python2的cupcake，应该克隆对应的**Py2\_v8.7.x.** 分支 首先得安装cogent环境 ```bash ## 创建anaCogent环境 conda create --name anaCogent python=2.7 anaconda ## 进入anaCogent环境 conda activate anaCogent ## 安装依赖 conda install -n anaCogent biopython -y conda install -n anaCogent -c bcbio bx-python -y conda install -n anaCogent -c conda-forge pulp -y ## 安装cogent cd git clone https://github.com/Magdoll/Cogent.git cd Cogent git checkout git submodule update --init --recursive cd Complete-Striped-Smith-Waterman-Library/src make export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/Cogent/Complete-Striped-Smith-Waterman-Library/src export PYTHONPATH=$PYTHONPATH:/Cogent/Complete-Striped-Smith-Waterman-Library/src cd ../../ python setup.py build python setup.py install ## 可以将export两行命令加入到.bashrc文件中 ## 开始安装cupcake，选择对应的Py2_v8.7.x. 分支 git clone -b Py2_v8.7.x https://github.com/Magdoll/cDNA_Cupcake.git cd cDNA_Cupcake/ python setup.py build python setup.py install ``` 6.2 安装依赖于python3.7环境的cupcake 参考 ```bash conda create --name Cupcake python=3.7 conda activate Cupcake conda install Biopython git clone https://github.com/Magdoll/cDNA_Cupcake.git ##添加到环境变量 export PATH=$PATH:/sequence/ export PATH=$PATH:/rarefaction/ ## 这不找不到，不加了 cd cDNA_Cupcake/ python setup.py build python setup.py install collapse_isoforms_by_sam.py -h ``` ### 6.3 具体的使用方法 参考 ```bash ## 将Gmap得到的sam结果文件进行排序 sort -k 3,3 -k 4,4n gmap.sam >gmap_sorted.sam ## 结合gmap的比对结果去除冗余的转录本 collapse_isoforms_by_sam.py --input polished.hq.fastq --fq -s gmap_sorted.sam -o tama/test --dun-merge-5-shorter -c 0.95 -i 0.85 ``` * \--input 输入文件 * \--fq 指定输入文件为fastq * -s Gmap输出后的sam文件经过sorted * -o 输出文件前缀，当然也可以加目录，直接输出到对应目录下 * -c 最小覆盖度 * -i 相似度 * \--dun-merge-5-shorter 5‘端的read由于测序的原因可能是真是存在差别，也可能是冗余；跟建库方式有关；因为设计引物数利用ployA的，所以5’端的序列可能没有完全扩到 **过滤因为5’端测序的误差，导致冗余没有完全去除** ```bash filter_away_subset.py test.collapsed ``` ## 7.与已有的注释信息进行比较 ```bash /usr/bin/python ~/software/MatchAnnot/matchAnnot.py --gtf 已经发表的基因组gtf文件 --format alt gmap的比对结果sam文件需要按照染色体顺序排好序 >111 ``` 报错 ```bash raise RuntimeError ('chromosome %s not in annotation' % chr) RuntimeError: chromosome scaffold_662 not in annotation ## 基因组注释文件里没有这个scaffold的煮熟信息，但是基因组序列在建库的时候有，比对的sam结果里就会有 ``` ## 8.Alternative splice.py脚本进行分类 ```bash git clone https://github.com/liangfan01/pipeline-for-isoseq.git --depth 1 ## 脚本在other目录下，依赖于python2 ## 先退出conda，回到bash下 conda deactive ## 安装依赖包 pip2 install svgwrite --user pip2 install networkx --user ## 我的python环境被污染了，只能指定对应的python /usr/bin/python alternative_splice.py -i gfffile -g 参考基因组文件gtf文件 -f 参考基因组文件 -o 输出路径 -os -as -ats T -op ``` :warning:alternative\_splice.py脚本中使用的参考基因组的gtf文件还需要使用awk，进行转化 ```bash awk '{print $1 "\t" $2 "\t" $3 "\t" $4 "\t" $5 "\t" $6 "\t" $7 "\t" $8 "\t" $11 " " $12 " " $9 " " $10}' cufflinks转化后的gtf文件 >最后可以使用的gtf文件 ``` ### 输出文件 ```bash $ tree . ├── acceptor.list.txt ├── alternative.splice.list.txt ├── donor.list.txt ├── error.orient.list.txt ├── gene.cluster.picture ├── novel.gene.list.txt ├── proved.transcript.list.txt ├── splice.ascode.list.txt ├── splice.ascode.stat.txt ├── transcript.cluster.list.txt ├── unproved.gene.list.txt └── unproved.transcript.list.txt ``` **暂时就更新到这里了\~\~\~** ## 鉴定可变剪切spladder 参考文档 在Ancona中安装这个软件 * 进行可变剪切的鉴定 ```bash spladder build -a ../../../07_annotation/merge.gtf -b ../test.bam -c 1 -o ./ --merge-strat single ``` ## 参考 完整的分析 流程 PacBio官方SMART软件使用说明**V8版本的** single molecular real-time GMap软件 samtools输出文件格式 全长转录本分类 Cupcake分析流程 去除冗余之后的分析流程 --- # Agent Instructions: Querying This Documentation If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question. Perform an HTTP GET request on the current page URL with the `ask` query parameter: ``` GET https://zpliu.gitbook.io/booknote/ke-bian-jian-qie/ruan-jian-shi-yong/01-san-dai-ce-xu-isoseq.md?ask= ``` The question should be specific, self-contained, and written in natural language. The response will contain a direct answer to the question and relevant excerpts and sources from the documentation. Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.