转录组ref流程工作手册.doc

转录组ref流程工作手册一、Reference 流程生物学原理1.1 实验流程图一：转录组实验流程当我们得到样品时，必须对其测序，才能得到分析所需的数据。测序基本过程：提取样品总RNA后，用带有Oligo(dT)的磁珠富集真核生物mRNA（若为原核生物，则用试剂盒去除rRNA后进入下一步）。加入fragmentation buffer将mRNA打断成短片段，以mRNA为模板，用六碱基随机引物（random hexamers）合成第一条cDNA链，然后加入缓冲液、dNTPs、RNase H 和DNA polymerase I合成第二条cDNA链，在经过QiaQuick PCR试剂盒纯化并加EB缓冲液洗脱之后做末端修复并连接测序接头，然后用琼脂糖凝胶电泳进行片段大小选择，最后进行PCR扩增，使用建好的测序文库进行测序。得到RNA的序列后，又可以找到它的参考序列（物种本身的基因、基因组）时，可以用reference流程对数据进行详细的分析。Reference后面所有的流程都是基于参考序列进行的，所以选择正确的参考序列十分重要。1.2信息分析流程 得到测序序列后，即可利用比对软件，将所测序列比对到参考基因或基因组上，并进行后续分析，信息分析流程图如下：图二：转录组信息流程1.2.1原始fq序列简介测序得到的原始图像数据经base calling转化为序列数据，我们称之为raw data或raw reads，结果以fastq文件格式存储，fastq文件为用户得到的最原始文件，里面存储reads的序列以及reads的测序质量。在fastq格式文件中每个read由四行描述：read IDTGGCGGAGGGATTTGAACCC+bbbbbbbbabbbbbbbbbbb每个序列共有4行，第1行和第3行是序列名称(有的fq文件为了节省存储空间会省略第三行“”后面的序列名称)，由测序仪产生；第2行是序列；第4行是序列的测序质量，每个字符对应第2行每个碱基，第四行每个字符对应的ASCII值减去64，即为该碱基的测序质量值，比如h 对应的ASCII值为104，那么其对应的碱基质量值是40。碱基质量值范围为0到40。表 1为Solexa测序错误率与测序质量值简明对应关系，具体计算公式如下：Qphred =-10 log10(e)表 1 Solexa测序错误率与测序质量值简明对应关系测序错误率测序质量值对应字符5%13M1%20T0.1%300.01%40h1.2.2原始fq序列处理 某些原始序列带有adaptor 序列，或含有少量低质量序列。我们首先经过一系列数据处理以去除杂质数据，得到Clean reads。按如下步骤进行处理：1. 去除含adaptor的reads2. 去除N的比例大于10%的reads3. 去除低质量reads（质量值Q *_final.sh.nohup &2.4 查看任务进展1） 操作任务的命令行：查看个人所有在跑任务： qstat -u * （用户是*）查看某一个在跑任务： qstat -j 24832|less （任务号是24832）杀掉个人所有在跑任务： qdel -u *杀掉某一个在跑任务： qdel 24832如果是因为某一个运行文件出错导致需要杀掉所要相关在跑任务，应该先杀掉这个在公共节点上跑的任务如：上面*_final.sh出错了，可以按以下步骤处理：top -u daichm按c键查看详细信息，找出所要杀掉的任务，假设*_final.sh对应的任务号是23849则可按k键，输入工作号，回车然后按9再回车即可杀掉该任务，再去做上面操作。2） 查看整个任务进展：a 查看*_final.sh.nohupb 进入part_shell目录，查看相应的任务运行信息，主要有可以查看以下几个文件：*. globle*.log进入下一层目录，查看.o和.e文件。找出问题所在并进行处理。2.5 任务完成1）结果检查：a,结题报告是否完整生成？ b,打包数据中，相关文件是否齐全？ c,分析要求是否都做好了，差异分析有没有漏掉？ d,有没有空文件产生？2）数据备份：由于各方面的原因，产生的数据有可能会丢失，建议对一些重要的数据在相对稳定的盘阵里做多一个备份，以免发生不必要的大麻烦。三、Reference 流程程序模块说明配置文件：ref.lib主程序脚本：maid.shperl reference_transcriptome_pipeline.pl -name huyang -lib ref.lib -outdir /ifs1/DGE_SR/daichm/project/HUYlfvT/result -diff -filter -2bwt -soap 2.20 -genome Populus_euphratica.0114.genome -gene Populus_euphratica.0114.cds -psl Populus_euphratica.0114.gff.psl -doall -verbose关键程序： reference_transcriptome_pipeline.pl其各项参数代表意思：Usage basic parameters:基本参数 -name species name(necessary!)即物种名，注意不是文库名 -lib input lib file(necessary!), a format example: file inputlib.配置文件 -outdir result dir(necessary!)结果输出目录 -genome genome sequence(necessary!)参考基因组 -gene gene sequence(necessary!)参考基因 -psl gene psl(necessary!)psl文件 analysis options:分析选项 -soap soap version(2.01|2.20|.)soap版本选择，目前用soap2.20 -filter Filter reads过滤数据，得到clean reads，一般也是必须选项 -div Divide analysis by chromosome name,if all chromosomes size is large.基因组大时按染色体分块处理 -doall Do all analysis below,including 5 parts.包括以下五个选项 -basic Do basic analysis.基本生物信息分析 -alter Alternative Splice analysis高级生物信息分析中的可变剪切 -novel Novel Transcript analysis高级生物信息分析中的发现新转录本 -utr ExtendGene analysis高级生物信息分析中的基因结构优化 -svg Produce SVG figure基因结构以及 Reads 在基因组上分布的精确图形 -diff Gene expression difference高级生物信息分析中的差异表达基因 -verbose output verbose information to screen输出运行信息到标准输出上 -help output help information to screen帮助文档分析的流程图：使用参数说明：命令示例：1) bwt,filter对应程序：/ifs1/DGE_SR/daichm/project/HUYlfvT/result/part_shell/bwt_filter.sh详细情况：a).基因组建库：/panfs/DGE_SR/hezengquan/soft/SOAPaligner/soap2.20release/2bwt-builder /ifs1/DGE_SR/daichm/project/HUYlfvT/Populus_euphratica.0114.genomeb).基因建库：/panfs/DGE_SR/hezengquan/soft/SOAPaligner/soap2.20release/2bwt-builder /ifs1/DGE_SR/daichm/project/HUYlfvT/Populus_euphratica.0114.cdsc).样本数据过滤(举其中一例)：sh /ifs1/DGE_SR/daichm/project/HUYlfvT/result/part_shell/Filter.huiyang_chuli_L1.sh2) soapa).对基因组所建的库跑soap：(举其中一例）/panfs/DGE_SR/hezengquan/soft/SOAPaligner/soap2.20release/soap -a /ifs1/DGE_SR/daichm/project/HUYlfvT/result/Reads/huiyang_chuli_L1_1.fq -b /ifs1/DGE_SR/daichm/project/HUYlfvT/result/Reads/huiyang_chuli_L1_2.fq -D /ifs1/DGE_SR/daichm/project/HUYlfvT/Populus_euphratica.0114.genome.index -m 0 -x 10000 -s 40 -l 35 -v 3 -o /ifs1/DGE_SR/daichm/project/HUYlfvT/result/Alignment/Genome/huiyang_chuli_L1.Genome.PESoap -2 /ifs1/DGE_SR/daichm/project/HUYlfvT/result/Alignment/Genome/huiyang_chuli_L1.Genome.PESoapSingleb).对基因所建的库跑soap：(举其中一例）/panfs/DGE_SR/hezengquan/soft/SOAPaligner/soap2.20release/soap -a /ifs1/DGE_SR/daichm/project/HUYlfvT/result/Reads/huiyang_chuli_L1_1.fq -b /ifs1/DGE_SR/daichm/project/HUYlfvT/result/Reads/huiyang_chuli_L1_2.fq -D /ifs1/DGE_SR/daichm/project/HUYlfvT/Populus_euphratica.0114.cds.index -m 0 -x 1000 -s 40 -l 35 -v 3 -r 2 -o /ifs1/DGE_SR/daichm/project/HUYlfvT/result/Alignment/Gene/huiyang_chuli_L1.Gene.PESoap -2 /ifs1/DGE_SR/daichm/project/HUYlfvT/result/Alignment/Gene/huiyang_chuli_L1.Gene.PESoapSingle关键程序：/panfs/DGE_SR/hezengquan/soft/SOAPaligner/soap2.20release/soap其各项参数代表意思：Usage: soap options -a query a file, *.fq, *.fa -b query b file -D reference sequences indexing table, *.index format -o output alignment file(txt) -M match mode for each read or the seed part of read, which shouldnt contain more than 2 mismaches, 4 0: exact match only 1: 1 mismatch match only 2: 2 mismatch match only 4: find the best hits -u output unmapped reads file -t output reads id instead reads name, none -l align the initial n bps as a seed 256 means whole length of read -n filter low-quality reads containing n Ns before alignment, 5 -r 0,1,2 how to report repeat hits, 0=none; 1=random one; 2=all, 1 -m minimal insert size allowed, 400 -x maximal insert size allowed, 600 -2 output file of unpaired alignment hits -v maximum number of mismatches allowed on a read. 5 bp -s minimal alignment length (for soft clip) 255 bp -g one continuous gap size allowed on a read. 0 bp -R for long insert size of pair end reads RF. none(means FR pair) -e will not allow gap exist inside n-bp edge of a read, default=5 -p number of processors to use, 1 -h this help3) posCoveragea)对单样本处理/ifs1/DGE_SR/daichm/project/HUYlfvT/result/part_shell/samples_pos.sh其中用到的程序为：/nas/DGE_SR01/daichm/ref/posCoverageb)合并所有样本/ifs1/DGE_SR/daichm/project/HUYlfvT/result/part_shell/posCoverage.sh其中用到的程序为：/nas/DGE_SR01/daichm/ref/merge_poscoverage.pl4) transcript-unit/ifs1/DGE_SR/daichm/project/HUYlfvT/result/part_shell/TranscritUnit.sha).PosCoverage.TAR其中用到的程序为：/nas/DGE_SR01/daichm/ref/Mask2Tar.plb).Filter其中用到的程序为：awk $335 print /ifs1/DGE_SR/daichm/project/HUYlfvT/result/Poscoverage/AllChr.AllTissue.PosCoverage.TAR /ifs1/DGE_SR/daichm/project/HUYlfvT/result/Poscoverage/AllTissue.PosCoverage.TAR.Filterc).PairEndJoinTAR其中用到的程序为： /nas/DGE_SR01/daichm/ref/PairEndJoinTAR.pld).TAR2Genes其中用到的程序为： /nas/DGE_SR01/daichm/ref/TAR2Genes.ple).TARGenes2psl其中用到的程序为：/nas/DGE_SR01/daichm/ref/TARGenes2psl.pl5) important analysis step/ifs1/DGE_SR/daichm/project/HUYlfvT/result/part_shell/analysis5.pla).BasicAnalysis&DiffBasicAnalysis:/ifs1/DGE_SR/daichm/project/HUYlfvT/result/part_shell/basic.sh其中用到的程序为：perl /nas/DGE_SR01/daichm/ref/MapReadsStat.pl perl /nas/DGE_SR01/daichm/ref/ReadsRandomInGene.plperl /nas/DGE_SR01/daichm/ref/Soap_Coverage.plperl /nas/DGE_SR01/daichm/ref/CoverageList.pl perl /nas/DGE_SR01/daichm/ref/chromosome_graph_wb.pl Diff:/ifs1/DGE_SR/daichm/project/HUYlfvT/result/part_shell/diff.sh其中用到的程序为：/nas/DGE_SR01/daichm/ref/DiffExp/DiffExp_pipeline.plb).Extend/ifs1/DGE_SR/daichm/project/HUYlfvT/result/part_shell/extend.sh其中用到的程序为：/nas/DGE_SR01/daichm/ref/getGene.pl/share/raid1/genome/bin/blat/nas/DGE_SR01/daichm/ref/pslbest.pl/nas/DGE_SR01/daichm/ref/TarBGFortholog.pl/nas/DGE_SR01/daichm/ref/UTR.plc).AlternativeSplicing/ifs1/DGE_SR/daichm/project/HUYlfvT/result/part_shell/GeneSpliceSite.sh其中用到的程序为：/nas/DGE_SR01/daichm/ref/GeneSpliceSite.pl/ifs1/DGE_SR/daichm/project/HUYlfvT/result/part_shell/JoinSplice.sh其中用到的程序为：/nas/DGE_SR01/daichm/ref/JoinTARForSoap.pl/panfs/DGE_SR/hezengquan/soft/SOAPaligner/soap2.20release/2bwt-builder/ifs1/DGE_SR/daichm/project/HUYlfvT/result/part_shell/TrimNomap.sh其中用到的程序为：/nas/DGE_SR01/daichm/ref/search_reads_TrimNomap.pl/ifs1/DGE_SR/daichm/project/HUYlfvT/result/part_shell/JunctionSoap.sh其中用到的程序为：/panfs/DGE_SR/hezengquan/soft/SOAPaligner/soap2.20release/soap/ifs1/DGE_SR/daichm/project/HUYlfvT/result/part_shell/AlternativeSplice.sh其中用到的程序为：/nas/DGE_SR01/daichm/ref/JunctionReadsStat.pl/nas/DGE_SR01/daichm/ref/SikppedEx