第四章Illumina測序數據分析方法簡介

第一節 下機數據的初步處理一、CASAVA的運行和參數設置在Illumina數據下機後的輸出文件中，並不直接存在後續分析所需要的Fastq文件，需要通過bcl2fastq Conversion Software（v1.8.4; support.illumina.com/ downloads/bcl2fastq_con version _software_184.html）（該軟體原名為CASAVA，後更名為bcl2fastq）來實現。（一）樣品清單

以處理Hiseq2500測序儀下機數據為例，在運行bcl2fastq時，需要提供所有樣本名稱與對應Index表格，稱為樣品清單（sample sheet），具體格式如圖4.1所展示；其中FCID為flow cell的編號；Lane列為樣本所在Lane編號；Sample ID列為樣本名稱；Index列為樣本對應的Index序列，同一條Lane的樣本不能使用相同Index；如果在同一條Lane中混合不同的測序樣本，建議盡量使每個樣本的Index與其他樣本之間至少有2個鹼基的差異；如果在HiSeq2500測序中引入了雙Index標記，在樣品清單中可以使用「-」連接兩組Index；Sample Project列為樣本所屬的研究項目名稱（project）。

圖4.圖4. * ARABIC1Hiseq 2500樣品清單（二）bcl2fastqbcl2fastq根據樣品清單表格中Sample Project構建一級目錄，同一Project下所有樣品放在下一級目錄。根據同一Project下樣本的Index不同，R1和R2端fastq文件分別輸出為壓縮格式，放在/Sample Project/Sample ID/目錄下。如果Index無法匹配則輸出到/Undetermined_ indice /Sample_lane/錄下。以處理Hiseq2500測序儀下機數據為例，bcl2fastq該命令的使用如下，包括兩步命令運行：perl/bcl2fastq-1.8.3/bin/configureBclToFastq.pl--input-dir $inputdir--output-dir$outputdir--sample-sheet $samplecsv--force--ignore-missing- bcl--use-bases-masky257n，I6n，Y257n --mismatches 1

make –j 10

參數中$input指向下機數據目錄；$ outputdir為輸出目錄；$samplecsv指向samplesheet；--mismatch 1參數為允許樣本Index有一個錯配。（三）其他測序儀的文件轉換2015年下半年，Illumina公司新推出了Hiseq3000/4000及Next500系列測序儀，由於測序原理上的改進，bcl2fastq Conversion Software2（v2.17）整體版本也進行了相應升級；樣品清單表格格式改為與Miseq系列的表格相似，具體情況和命令變化見圖4.2。

圖4.圖4. * ARABIC2Hiseq 3000/4000樣品清單相應的運行命令更改為：/bcl2fastq2-2.17/bin/bcl2fastq -i $inputdir --output-dir $outputdir --sample-sheet $samplecsv--ignore-missing-bcl --use-bases-mask y150n，I6n，Y150n --ignore-missing-filter

--ignore-missing-controls --ignore-missing-positions

（四）文件轉化前後的目錄結構圖4.3中所展示的為下機數據原始目錄結構和通過bcl2FASTQ處理後的文件目錄。

圖4.圖4. * ARABIC3原始數據和處理後數據結構對比（五）查看下機數據產量和質量（1）通過bcl2fastq生成的文件查看在完成bcl2fastq後可以通過生成目錄下文件/Basecall_Stats_[FCID]/Demultiplex_ Stats.htm文件查看本次數據的整體產量和質量，如下表。

（2）通過SAV查看器查看在運行bcl2fastq之前，同樣可以獲得本次數據的產量與質量，通過Illumina Sequencing Analysis Viewer（SAV，support.illumina.com/se _analysis_viewer_sav.html）獲得；通過SAV查看器打開含有InterOp、RunInfo.xml、runParame ters.xml 3個文件的文件夾（這3個文件位於原始下機目錄中），如圖4.4、圖4.5所示。

圖4.4 SAV文件展示1

圖4.5 SAV文件展示2二、下機數據的質控處理

（一）Fastq文件格式

由於實驗操作和測序儀器等原因，會導致測序數據中部分短序列（reads）尾部質量下降或接頭（adapter）自連，這些序列會對後期的數據分析造成困擾，所以在拿到原始的Fastq文件後我們需要對數據進行質量控制，去掉一些低質量的序列。首先我們需要了解下機數據的Fastq文件格式；Fastq文件為4行一組的序列格式。@SEQ_IDGATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT+!*（（（（***+））%%%++）（%%%%）.1***-+*））**55CCF>>>>>>CCCCCCC65@列為序列名字，其中包含一個絕對名字、樣本所在Lane、Cluster坐標位置和R1/R2端標誌等信息。例如，@HWI-D00447:78:XXXXXXX:1:1101:1217:1986 1:N:0:ACAGTG；其中HWI-D00447:78: XXXXXXX為絕對名字；1表示位於Lane1；1101為晶元號；1217:1986為坐標位置；1:N:0: ACAGTG表示該條序列為R1端序列，N表示序列沒有經過CASVA filter（如果設置過CASAVA過濾reads，這裡顯示Y）；ACAGT表示該條序列的Index。第二行為具體序列信息。第三行為「+」。

第四行為該條序列每個鹼基的質量值（計算公式如下），表示為該鹼基的錯誤率：將轉換後的值通過ASCII轉換為字元後顯示；這種方式被稱為Phred；現在HiSeq測序儀使用phred+64的計算方式，數據質量從0到40；分別對應!"#$%&（）*+，-./0123456789:;<=>?@ABCDEFGHI；質量值為10時說明，該鹼基的錯誤率為10%；20時表示錯誤率為1%；40時表示錯誤率為0.01%。

（二）通過Trimmomatic進行原始數據質控根據每條序列中帶有的質量信息，我們可以通過一定的過濾條件去掉一些質量過低的序列和鹼基；如使用Trimmomatic（Bolger et al.，2014）（usadellab.org/cms/?）、fastx（Gordon and Hannon，2010）（http:// hannonlab.cshl.edu/fast）等軟體實現。Trimmomatic可以非常方便地實現一次運行多個參數的過濾條件，運行參數如下：java -jar trimmomatic-0.30.jar PE -phred33 input_forward.fq.gzinput_reverse.fq.gzoutput_forward_paired.fq.gz output_forward_unpaired.fq.gzoutput_reverse_paired.fq.gz

output_reverse_unpaired.fq.gz

ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36命令中通過給定雙端測序的原始fastq.gz文件input_forward.fq.gz和input_ reverse.fq.gz，指定輸出文件為4個文件，分別為R1配對和未配對序列文件，R2端配對和未配對序列文件；同時設定接頭（adapter）序列文件TruSeq3-PE.fa（軟體自帶，該文件為Illumina連接序列用的固定序列）；並制定了以下數據過濾原則：如序列頭部和尾部質量值都低於3的鹼基去掉，以4個鹼基為窗口滑動，將平均質量低於15的窗口以後的鹼基丟掉，最後設定保留最小長度為36 bp。（三）FastxFastx軟體更加靈活，可以通過15個子程序分別實現序列的切割和轉換，根據低質量比例去除reads；如去除質量值低於20的鹼基占序列總長度50%的序列。fastq_quality_filter -q 20 -p 50 –i R1.fastq -o R1_filter.fastq（四）對比原始序列和過濾後序列情況圖4.6為對比原始序列和過濾後序列情況，可以通過Fastqc（bioin

圖4.6 fastqc鹼基質量展示

formatics. babraham.ac.uk/projects）軟體對比序列質量。圖中橫坐標表示鹼基位置，縱坐標表示每個位置上鹼基質量的情況，左圖表示質控前，右圖表示質控後。（高 強）第二節DNA測序數據分析簡介一、基因組從頭測序數據分析基因組從頭（de novo）測序是指在沒有參考基因組的情況下，對物種的基因組進行測序、拼接和組裝，進而獲得物種的基因組序列圖譜。物種基因組序列圖譜的完成，是進一步研究該物種的遺傳信息與進化的基礎。2010年熊貓基因組（Li et al.，2010）的完成，標誌著二代從頭測序技術在國內的開始，採用全基因組鳥槍法（whole genome shotgun，WGS）與二代測序相結合的技術，構建了BAC文庫和插入片段長度分別為500 bp、2 kb、5 kb、10 kb的二代測序文庫，獲得平均讀長52 bp的序列，可用鹼基約176 GB（測序深度約為73×），最終組裝出了大熊貓的21條染色體，約2.4 Gb的基因組序列。contig N50長約40 kb，同時注釋了2萬多個基因。通過同源基因的分析發現大熊貓不喜歡吃肉的主要原因是T1R1基因失活，無法感覺到肉的鮮味。類似的，2013年發表的藏豬基因組同樣也是利用二代測序技術完成了基因組從頭測序。相對而言，某些植物和昆蟲的基因組龐大、複雜度高、重複序列多，從頭測序難度更大，但是時隨著二代測序技術的提升，如序列讀長的增加，測序時間的降低，一些基因組複雜的多倍體物種，如二倍體小麥A基因組（Ling et al.，2013）和D基因組草圖（Jia et al.，2013）已經完成。伴隨著讀長更長的新一代測序技術的發展，相信未來會有越來越多的物種基因組得到解析，進而讓我們更加深入地了解物種進化的歷史。（一）分析流程基因組從頭測序序列的組裝是項複雜的工作，在計算過程中涉及大量的數據和計算資源消耗，在此我們僅列出一個基本的分析流程（圖4.7）和簡要的分析參數。在測序數據方面，不同複雜度的基因組數據測序深度需100X～200X。（二）分析參數

要進行基因組從頭（de novo）測序，研究者首先要對物種進行基因組調查（survey），了解基因組的大小和雜合度等情況。基因組大小可以通過查詢相關資料庫獲得，如genomesize.com/；或者通過實驗獲得，如流式細胞儀（Yoshida et al.，2010）測定，福爾根染色，定量PCR和k-mer估計等。

圖4.7從頭測序分析流程估算k-mer值不僅可以估算基因組大小，同時也能檢測基因組的雜合度，雜合度越高的物種基因組組裝難度越大。做k-mer估計時常用軟體有Jellyfish（Lynam et al.，2006）（www. cbcb. umd.edu/ software /jellyfish/）和KmerGenie（Chikhi and Medvedev，2013），通過分析獲得k-mer的統計頻數，如圖4.8所示。

圖4. 8 k-mer估計圖4.9為牡蠣基因組（Zhang et al.，2012）在17-mer中調查發現的情況。雙峰的情況說明基因組雜合度很大，超過35.49%的17-mer測序深度超過255×，說明基因組序列中存在大量的重複序列。

圖4. 9基因組雜合度可以通過k-mer的峰值情況（三）SOAPdenovo目前有很多基因組組裝軟體工具可供選擇，如DNA序列組裝常用的SOAPdenovo（Xie et al.，2014a）（soap.genomics.org.cn/so）、MaSuRCA（Zimin et al.，2013）（gen ome .umd.edu /masurca.html）等軟體。SOAPdenovo的優點是組裝速度非常快，對於內存和磁碟空間消耗較小，但缺點是所得的組裝序列較短。MaSuRCA是2013年發表的一種功能強大的軟體工具，軟體自身帶有鹼基質量等處理過程，得到的contig序列長，但是組裝耗費時間長，組裝過程對磁碟空間消耗也大。如果我們用上述兩種軟體分別對10 G水稻數據進行組裝，SOAPdenovo在15～30 min內即可結束，得到的contig N50長度約400 bp；而MaSuRCA耗時約為24 h，contig N50約為4 kb。SOAPdenovo組裝運行命令：SOAPdenovo all -p 8 -s config.file -K 31 -R -o Sample 1 > ./Sample.log 2> Sample.errConfig.file中包含參數意義如下：avg_ins=280：表示平均插入片段reverse_seq=0asm_flags=1：組裝步驟1表示contig；2表示scaffoldrank=1q1=R1.fastqq2=R2.fastq我們一般通過N50/N90等指標來評價組裝結果，具體可以通過abyss（Simpson et al.，2009b）帶的一個abyss_fac.pl腳本進行，具體命令如下。perl abyss_fac.pl Sample.ctg.fasta >> Sample.ctg.result二、重測序數據分析流程對於完成從頭測序的物種，我們已經獲得了基因組的完整序列信息，利用全基因組重測序技術對其個體或群體的基因組進行測序及差異分析，可獲得SNP、InDel、SV、CNV等大量的遺傳變異信息，進而可以對該物種的基因功能挖掘和群體進化進行深入分析。隨著2002年人類基因組測序工作的完成，人類基因組重測序已經成為人類遺傳學和轉化醫學的重要研究手段。在2012年發表的自閉症（Jiang et al.，2013）研究中，通過對32個自閉症家系進行全基因組重測序（30X深度），對單鹼基突變和拷貝數變化進行檢測，在綜合分析已經發現的9個自閉症相關基因的基礎上，進一步發現了4個新的自閉症基因和8個候選自閉症風險基因。近兩年隨著測序成本的下降，重測序技術已經深入到人們的日常生活中，如對於胎兒的21三體綜合征、自閉症等孕期篩查和通過marker對癌症進行早期篩查。在科研領域，對重要的農作物進行重測序，對於發覺其基因功能也是至關重要的。2014年發表的水稻（Chen et al.，2014a）研究中，選取529份水稻材料進行全基因組重測序，對兩個亞群水稻代謝性狀進行全基因組關聯分析（GWAS），鑒定出2947個與634個基因相關的主導SNP位點。隨後，在210個秈稻的自交群體中進行驗證，定位出36個候選基因與代謝相關。對36個候選基因進行實驗驗證，最終確定了5個候選基因。2015年發表的302份代表性大豆材料（Zhou et al.，2015）的全基因組重測序工作表明：大豆在馴化過程中受到了強烈的選擇瓶頸效應，鑒定出121個強選擇信號；同時對大豆的種子大小、種皮顏色、生長習性、油含量等性狀做了全基因組關聯分析，找出了一系列顯著關聯位點；把選擇信號、GWAS信號及前人研究的油含量QTL相整合，發現很多選擇信號和油相關性狀有關，說明大豆產油性狀受人工選擇較多，形成了複雜的網路系統，共同調控油的代謝，從而引起不同種油質相關性狀的變異。研究還定位了一些重要農藝性狀的調控位點，並且明確了一些基因在區域化選擇中的作用，如控制花周期的E1、控制生長習性的Dt1、控制絨毛顏色的T等。這為大豆重要農藝性狀調控網路的研究奠定了重要基礎。除了利用全基因組鳥槍法（WGS）對基因組進行隨機打斷測序，獲得全基因組的序列信息外，近年來還出現了成本更加低廉的重測序方法，如GBS、RAD等簡化基因組方法，通過對特定酶酶切後的基因組DNA片段進行高通量測序，可以保證在更低數據量的情況下獲得全基因組的SNP信息。例如，2013年發表的針對於高粱株高（Morris et al.，2013）的研究中，通過對971份高粱材料進行GBS高通量測序，產生了21 G的測序數據，獲得265 487個SNP；對其中336株高粱個體進行全基因組關聯分析，發現多個與株高相關的已知基因，並定位到多個與花序結構相關的候選基因。（一）重測序分析流程基因組重測序分析的一般流程如圖4.10所示

圖4.10重測序分析流程（二）mapping群體分析的第一步是將測序得到的短序列比對（mapping）回基因組，常用的DNA比對軟體有BWA（Li and Durbin，2009b）。1）BWA比對之前需要對參考基因組構建索引（index）。bwa inde reference.fasta表示對參考基因組構建索引。bwa mem –t 8 ref.fa reads.fq > aln-se.sam表示單端（single end，SE）測序數據，通過8個核運行。bwa mem –t 8 –M ref.fa read1.fq read2.fq > aln-pe.sam表示雙端（pair end，PE）測序數據，通過8個核運行，對短序列進行處理。2）文件格式轉換和變異檢測：比對生成的文件為sam文件，可以通過Samtools（Li et al.，2009）進行進一步處理和查看，其中包含了每條序列的比對位置、比對的具體情況；sam格式一般較大，通過轉換一般存儲為二進位壓縮的bam格式；並且比對後的序列，根據比對位置進行排序後，可以進行可視化查看和變異檢測。具體可以通過如下命令進行。smatools view -@ 8 –bS aln-pe.sam >> aln-pe.bam表示sam壓縮為bam；-@為多核參數。samtools sort -@ aln-pe.bam aln-pe.sorted.bam表示bam文件排序。（三）重測序的數據質量評估完成比對得到bam文件後，通過分析bam文件我們可以對重測序的數據質量進行評估，常見的評估指標包括：①比對短序列比例（mapping ratio），即比對上的短序列佔總測序序列的比例；②深度（depth），即基因組上被覆蓋的鹼基平均被多少條短序列（reads）覆蓋；③覆蓋度（coverage），即基因組多少鹼基被reads覆蓋；④插入片段長度（insert size），即指建庫時打斷片段長度，通過R1和R2端序列比對在參考基因組上的間距，構建頻數分布圖，一般檢測峰值長度和是否是單峰。1）運行samtools flagstat aln-pe.sorted.bam命令，屏幕輸出短序列比例，結果示例如下。

2）運行soap.coverage -cvg -sam -p 5 -i pe.sam -refsingle ref.fasta –o pe.coveage命令，屏幕輸出覆蓋度和測序深度，結果示例如下。

3）運行java-jar picard/CollectInsertSizeMetrics.jar INPUT=aln-pe.sort.bam DEVIATIONS=10 HISTOGRAM_FILE= pe.hist.DV10.pdf OUTPUT= pe.DV10.txt命令，輸出文庫不同插入片段長度頻數分布圖，示例結果如圖4.11所示。

圖4.11插入片段長度頻數分布圖同時，我們也可以通過IGV等可視化軟體，查看短序列在參考基因組的比對情況。（四）變異檢測對於變異檢測常見的是尋找個體間的單核苷酸多態性（SNP）和一些短的插入缺失突變（InDel）；常用的軟體包括GATK（McKenna et al.，2010a）、Samtools（Li et al.，2009）和SOAPsnp（Li et al.，2008）。1. GATK（1）GATK分析流程基因組分析工具包GATK（genome analysis toolkit，www. broadinstitute. org/gatk/）是由BROAD中心開發的一套專門針對高通量測序數據進行變異檢測的工具，該軟體在人和動物群體中應用相當廣泛，該軟體提供了詳細的操作規範和使用流程信息（broadinstitute.org/gatk data-processing -ovw），具體分析流程如圖4.12所示。（2）GATK運行命令GATK軟體推薦使用Bwa mem進行短序列比對，通過Picard（broadinstitute.github.io）進行重複序列的標記，然後根據流程調用GATK工具的不同方法進行操作，具體流程可參考以下腳本。

圖4.12 GATK推薦分析流程1）#輸入排序後的Bam文件。java -Djava.io.tmpdir=./tmp -jar picard/MarkDuplicates.jar I=PE.sorted.bamO=/PE.marked.bam METRICS_FILE=/PE.metricsFile CREATE_INDEX=trueVALIDATION_STRINGENCY=LENIENT REMOVE_DUPLICATES=trueASSUME_SORTED=truejava -Djava.io.tmpdir=./tmp -jar GenomeAnalysisTK.jar -nt 4 -T RealignerTargetCreator-R reference.fasta -o ./PE.intervals.list -I ./PE.marked.bamjava -Djava.io.tmpdir=./tmp -jar GenomeAnalysisTK.jar -T IndelRealigner –Rreference.fasta -I ./PE.marked.bam --targetIntervals ./PE.intervals.list -o ./PE.realigned.bamjava -Xmx20g -Djava.io.tmpdir=./tmp -jar GenomeAnalysisTK.jar -nct 4 –TBaseRecalibrator -l INFO -R reference.fasta -I ./PE.realigned.bam –covReadGroupCovariate -cov QualityScoreCovariate -cov CycleCovariate -o ./PE.b.recalFile--run_without_dbsnp_potentially_ruining_quality2）#單個個體的變異檢測。java -Djava.io.tmpdir=./tmp -jar GenomeAnalysisTK.jar -nct 4 -T PrintReads –Rreference.fasta -I ./PE.realigned.bam -o ./PE.recalibrated.bam -BQSR ./PE.b. recalFilejava -Xmx15g -Djava.io.tmpdir=./tmp -jar GenomeAnalysisTK.jar -T HaplotypeCaller –Rreference.fasta -I PE.recalibrated.bam -o ./raw.vcf --genotyping_mode DISCOVERY3）對於單個個體的變異檢測可以根據一步直接獲得；而對於群體數據需要先分別對個體樣本進行以上腳本。a）#對Sample1獲得GVCF文件。java -Djava.io.tmpdir=./tmp0 -jar GenomeAnalysisTK.jar -T HaplotypeCaller -R reference.fasta -I Sample1.recalibrated.bam -o ./Sample1.gvcf --genotyping_mode DISCOVERY --emitRefConfidence GVCFb）#對Sample2獲得GVCF文件。java -Djava.io.tmpdir=./tmp0 -jar GenomeAnalysisTK.jar -T HaplotypeCaller -R reference.fasta -I Sample2.recalibrated.bam -o ./Sample2.gvcf --genotyping_mode DISCOVERY --emitRefConfidence GVCFc）#合併GVCF文件得到群體變異數據。java -Djava.io.tmpdir=./tmp -jar GenomeAnalysisTK.jar -T GenotypeGVCFs -R reference.fasta -V Sample1.gvcf -V Sample2.gvcf -o ./combine.vcf -nt 8對於目前通用的變異文件輸出文件，大部分軟體都兼容VCF格式，目前已經更新到4.2版本（samtools.github.io/hts-）；具體格式如下。

VCF文件主要包括頭部「#」開始的注釋文件信息，包括格式版本號、使用的參考基因組地址和參考基因組染色體長度等；從「#Chrom」開始，為VCF文件正文，文件前9列為固定格式，分別為如下信息：檢測到變異的染色體、位置、ID、參考基因組鹼基信息、觀察到個體的鹼基信息、質量值、過濾標籤、信息列（包括總深度，觀察到的樣本等該位點的具體信息）和樣本標籤格式；從第十列開始為每個樣本的具體情況。2. SamtoolsSamtools同樣也能快速地檢測變異，並且使用方便，同時生成兼容的VCF格式。基本操作如下。#通過Samtools調用mpileup的方法，對Sample1和Sample2進行群體變異檢測，同時通過bcftools轉換為vcf格式samtools1.2 mpileup -t DP，DP4 -ugf reference.fasta Sample1.sorted.bam Sample2.sorted.bam |bcftools1.2 call -vmO z -o samtools.vcf.gz群體測序的關聯分析過程三、群體分析流程在群體分析中，通過高通量測序獲得群體基因組變異數據後，我們可以與表型數據結合進行基因定位工作，也可以通過野生品種與現有品種比較進行馴化與進化分析。在獲得群體變異數據後一般可以進行以下分析（圖4.12）。

圖4.12群體分析流程（二）SNP注釋在獲得群體基因組變異情況後，可以通過ANOVA、SNPEFF等軟體對VCF文件進行注釋，分析變異是否發生在基因區域，深入分析是同義突變還是非同義突變。在植物研究中，通過輻射手段進行誘變育種，可以快速獲得不同性狀個體。例如，對水稻種質資源Oryza sativaL. ssp.indica（9311）（Belfield et al.，2012）輻射誘變突變體Red-1進行20倍高通量測序，發現9.19%的基因組序列發生了突變，其中包含381 403個SNP、50 116個長度1～5 bp的InDel和1279個拷貝數變異，共涉及14 493個基因。在獲得群體或者個體變異信息後，我們可以通過以下命令獲得變異的注釋信息：java -Xmx4g -jar snpEff.jar eff -c snpEff.config -v rice_database filter.vcf > anno.vcf（三）基因型與表型關聯分析在獲得群體基因組變異數據後，可以結合個體表型數據（如身高、體重、患病、株高、粒重等）進行連鎖不平衡關係分析，獲得與表型性狀相關的候選基因區域，這種方法稱為全基因組關聯分析（GWAS）。利用GWAS技術，2005年首次進行了人類年齡相關性（Klein et al.，2005）黃斑變性的研究。後續通過這種方法在人類疾病方面發現了很多致病基因，如自閉症基因、囊性纖維化基因、亨廷頓病基因等。在植物方面，GWAS方法也定位了很多與植物重要性狀（如產量、株高等）相關的候選基因，在很大程度上加速了育種的進程。在分析方面很多軟體可以實現表型與基因型聯合分析，如人類基因組方面常用的plink（Purcell et al.，2007）（pngu.mgh.harvard.edu/～purcell/plink/），植物中常用的Tassel（sourceforge. net/projects/tassel/）和EMMAX（Zhou and Stephens，2012）（genetics.cs.ucla.edu/em）等軟體。（1）格式轉換在分析中，首先要進行文件格式的轉化，大部分的變異檢測結果都以VCF文件格式存在，首先通過VCFtools（Danecek et al.，2011）將VCF文件轉換為ped格式或者TPED格式，具體命令如下。vcftools --vcf all_last.vcf –plink #轉換為plink PED格式vcftools --vcf all_last.vcf --tped-plink #轉換為plink TPED格式（2）SNP過濾對於轉換後的文件一般進行SNP過濾，基本過濾條件為：MAF（最小等位頻率）>0.05，Gene>0.05（對於單個位點檢測到的樣本數目大於95%），mind>0.05（個體基因型數據缺失不超過5%）。plink --file out --noweb --maf 0.01 --geno 0.05 --mind 0.05 --out gwas #過濾後以ped格式輸出（3）關聯分析再與表型數據聯合分析得到每個SNP位點與基因型的關聯程度，操作命令如下：plink --file gwas --1 --pheno pheno.txt --all-pheno --assoc –adjust結果文件通過p值（顯著值）表示該位點與該表型的相關程度，一般我們通過繪製Manhattan plot（曼哈頓圖）來展示結果（圖4.13）；圖中橫坐標為12條染色體，縱坐標為-log10（p）。（四）QTL定位在目標性狀定位中，2012年發表的MutMap（Abe et al.，2012）方法，以野生型親本為參考，通過全基因組比對，計算SNP頻率，篩選SNP-index=1的位點，找到與表型變化相關的SNP區域。通過野生型親本和子代DNA池SNP頻率差異

染色體圖4.13 GWAS結果曼哈頓圖展示定位分析，找到葉片淺綠突變、莖長、葉長、花序數目、花序長度等相關的突變位點。2013年發表了一篇通過混池分組分析法（bulk segregant analysis，BSA）定位極端性狀（Takagi et al.，2013）的方法，將20～50株通過EMS誘變產生的極端性狀子代個體分別混合成2個DNA池，對其及親本分別建庫，利用Illumina平台進行基因組重測序，通過檢測全基因組的變異，分析頻率差異，成功定位了SNP-index差異群體重要農業性狀的QTL。目前這種方法在植物RILs群體和Fx群體中得到了廣泛使用。該軟體作者也發表了相關軟體和操作流程，可參考QTL-seq_framework1.4.4和MutMap_framework1.4.4（genome-e.ibrc.or.jp/ home/bioinformatics-team/mutmap/）。QTL-seq的結果如圖4.14所示。

圖4.14 QTL-Seq結果展示（另見彩圖）圖4.14為QTL-seq的展示結果，圖中點表示高池與低池的SNP頻率差異，頻率差異大的地方可能為極端表型的控制區域，通過滑動窗口的方法平均獲取一定區域的差異值，通過紅線表示出來，可以看到紅線在7號染色體中部表現出極端情況，所以預測該區域與表型相關。高通量測序同樣也適用於傳統的QTL定位。2013年，Huang等通過BinMap（Huang et al.，2009）方法構建了高密度遺傳圖譜進行QTL定位：通過對150個水稻RIL群體進行全基因組重測序，得到150萬個SNP，利用滑動窗口來確定染色體各個區段的歸屬，通過群體的交換位點，構建出高密度遺傳圖譜進行QTL定位（圖4.15）。這種方法大大提高了定位的精度和準確度，大大縮短了育種周期。此外，BinMap方法也成功應用於710份玉米F2群體（Chen et al.，2014b）GBS測序研究中，利用檢測到的1155158個SNP，構建了一個含有6533個bin-maker

圖4.15 BinMap結果展示（另見彩圖）的遺傳圖譜，遺傳圖距長度為1396 cM。通過對2個已知基因的準確定位，驗證了該圖譜的高質量和高準確性。對控制群體雄穗分枝數、穗行數及雌穗長度的區域進行定位，得到10個QTL，其中7個與前人報道的QTL相重疊，有3個MADS- box蛋白和一個BTB/POZ蛋白編碼基因位於qTBN5和qTBN7（分別長800 kb與1.6 Mb）之中，可能參與了雄穗結構的形成。（高 強）第三節 轉錄組測序標準信息分析轉錄組測序（RNA-seq）分析根據其是否依賴參考基因組信息，主要分為兩大部分：轉錄組從頭測序（de novosequencing）測序及重測序（re-sequencing）。轉錄組從頭測序是在不依賴物種基因組序列信息的前提下，用新一代高通量測序技術對某物種的特定組織或者細胞的轉錄本進行測序並獲取相應的轉錄本信息的過程。轉錄組重測序針對的是具有參考基因組序列的物種，用新一代高通量測序技術對某物種的特定組織或者細胞進行轉錄組測序，並與參考基因組進行比較，從中得到基因表達差異、可變剪接、融合基因等轉錄調控信息的過程。一、轉錄組數據分析（一）數據質控在進行轉錄組的從頭測序和重測序時，首先要做的是數據質控處理，即對原始數據進行去除接頭污染序列及低質量序列（reads）的處理。轉錄組測序數據與DNA測序數據的質控相似，在數據從Illumina測序儀下機之後，利用Illumina提供的軟體包bcl2fastq Conversion Software進行圖像數字信息轉換，產出Fastq格式數據。再通過Trimmomatic、CutAdapter及Fastx（hannonlab.cshl.edu/fast toolkit/）軟體包去除接頭污染序列，並對其低質量鹼基進行掃描過濾，從而得到有效數據（clean data）。這一部分序列處理與DNA測序數據處理過程是相同的，具體處理過程請詳見本章第一節「下機數據的初步處理」。（二）轉錄組序列組裝及比對經過質控之後，RNA-seq數據分析的第一步就是要把那些測序得到的短序列（reads）比對到該物種參考基因組或者轉錄組序列上。這取決於被測序物種是否具有高質量的參考基因組序列。如果沒有參考基因組，則需要將所得測序序列組裝成轉錄本序列後，再與之進行比對。轉錄組測序可用於未知基因組物種。通過軟體組裝可以獲得相應的轉錄本序列信息，用於研究基因結構、基因功能、可變剪接和新基因預測等。值得注意的是，轉錄組測序數據針對的是基因組中特定的基因區域，因此對於等量的測序數據來說，轉錄組測序數據的基因組覆蓋率要遠高於DNA測序數據。目前對於測序數據量，許多測序服務商都給出了各自的推薦標準。一般來說，更多的測序數據可以更好地保證拼接組裝的完整性。RNA-seq數據的組裝類似於基因組數據的組裝。短讀長的基因組裝配工具，如Velvet（Zerbino，2010）、ABySS（Simpson et al.，2009a）和ALLPATHS（Butler et al.，2008），不能直接應用於轉錄組組裝。原因如下：①在基因組測序中，DNA測序深度的預計值是基本一致的，而轉錄本的測序深度卻常常相差幾個數量級；②由於可變剪接是比基因組組裝中的線性問題更複雜的一個轉錄組組裝問題，通常需要作圖來表示每個位點的多個可變轉錄本。這些特點使得轉錄組裝配比基因組組裝的計算問題更加複雜多變。用於組裝轉錄本序列的軟體有很多，包括SOAPdenovo-Trans（Xie et al.，2014b）、Trinity（Haas et al.，2013）、Trans-ABySS（Robertson et al.，2010）等。這些軟體各具特點，逐一簡介如下。1. SOAPdenovo-TransSOAPdenovo-Trans由華大基因開發，用於組裝轉錄組數據，最新的版本下載地址為：sourceforge.net/project，版本為1.03，下載包中包括兩個k-mer版本的程序，即「SOAPdenovo-Trans-31mer」和「SOAPdenovo-Trans-127mer」。k-mer是序列組裝軟體中一個非常重要的參數，定義兩個reads之間重疊的長度，用于衡量兩個reads是否是連續的序列。k-mer值越大，就表示高表達的轉錄本能夠更完整地組裝出來（Surget-Groba and Montoya-Burgos，2010），與此相反，低表達的基因有可能在較小的k-mer值條件下組裝得更加準確。因此，k-mer值的選取取決於轉錄組拼接的具體需求：如果要得到更多樣化的轉錄本，應該適當降低k-mer值；如果要得到更長的轉錄本，就需要提高k-mer值進行組裝。一般來說，選取一個中間的k-mer值可以用來平衡兩個極端的組裝效果。使用SOAPdenovo-Trans軟體進行組裝的命令通常如下。$SOAP_HOME/SOAPdenovo-Trans-127mer all -s $config_file -o prefix_name -K 51 -p 40 -d 2 -e 5 -M 1 -L 200 -t 5 -G 5表明使用SOAPdenovo-Trans-127mer在k-mer長度為51的條件下，使用40個CPU進行組裝，「-d」命令通過k-mer的頻率去除錯誤reads，而「-e」、「-M」命令均為構建de Bruijin圖中的參數。「-L」、「-t」及「-G」命令用於組裝後的contig和scaffold的長度、數目及gap長度的過濾。Config_file中主要包含了組裝參數及數據路徑等信息，如圖4.16所示。

圖4.16 SOAPdenovo-Trans的參數配置文件SOAPdenovo-Trans軟體的優點是速度快，內存消耗較小。2. TrinityTrinity是專門為轉錄組的組裝設計的一種工具（Haas et al.，2013）。它首先將單個RNA-seq讀長擴展至更長的contig，然後用這些contig構建許多de Bruijn圖，然後在每幅圖中得到所有的剪接異構體代表路徑。命令非常簡單，即：$TRINITY_HOME/Trinity.pl --seqType fq --single single.fq --JM 20G$TRINITY_HOME/Trinity.pl --seqType fq --left left.fq --right right.fq --JM 20G對於鏈特異性的RNA-seq數據來講，需要加上額外的參數「--SS_lib_type」用來區分不同的建庫方法，如圖4.17所示。

圖4.17 Trinity對於不同建庫方法的參數選擇組裝後得到的轉錄本序列可以通過Trinity自帶軟體TrinityStats.pl進行信息統計，命令如下。$TRINITY_HOME/util/TrinityStats.pl trinity_out_dir/Trinity.fasta運行後可以得到轉錄本及基因的數目、N50、contig的平均長度等信息。N50是評價基因組組裝質量的一個常見指標，即組裝序列按長度從大到小進行排序後，取覆蓋總長度50%時的contig長度。這個值越大，表明組裝的序列整體長度越長，組裝的效果越好。而在轉錄組組裝過程中，Ｎ50不太適用。評價轉錄組組裝質量時，可以確定一套已知參考基因序列集合（同一物種或相近物種），將組裝序列與其進行比較，估計組裝序列在參考基因集合中的覆蓋率，以及覆蓋全長參考基因集合的數目。這樣得到的統計數值更有意義。Surget-Groba和Montoya-Burgos（2010）提出了scaffolding using translation mapping（STM）的組裝方法，通過翻譯contig序列並與參考蛋白質序列進行比對，將比對到相同參考蛋白質序列的contig進行匯總並進一步組裝，提高組裝的正確性，從而提高組裝質量。對於組裝完成的序列，可以進一步通過長度、氨基酸讀碼框進行篩選。篩選得到的具有較長可讀框（ORF）的序列稱為Unigene集合，可以用於進一步的分析。（三）基因表達分析及差異表達基因的篩選轉錄組分析中最基礎和最常規的工作是基因差異表達分析。分析時，首先要確定基因的表達量，再進行樣品間或基因間的表達量差異分析。獲取基因的表達量，首先要將基因比對到參考基因組序列上。如果被測序物種已經具有較高質量的參考基因組序列，就可以將reads比對到參考基因組序列上。與DNA測序數據不同，轉錄組測序數據比對到參考基因組上時，需要考慮參考基因組中內含子（intron）區域對比對過程的影響。因為轉錄組中不含有內含子，所以比對軟體需要能夠處理大片段的缺失（gap），用以跨過參考基因組中的內含子區域。目前常用的轉錄組比對軟體有Tophat（Trapnell et al.，2009）、Star（Dobin et al.，2013）、Hisat（Kim et al.，2015）等。對於未知基因組的物種來說，可以把組裝後的轉錄本序列作為參考序列來進行比對，比對軟體除了同樣採用Hisat等軟體之外，由於組裝的轉錄本序列不含有內含子，因此，也可以採用DNA序列比對軟體，如BWA（Li and Durbin，2009a）等。Hisat在運行時要構建大量的片段索引，不僅可以減少內存消耗，而且其比對速度要遠快於Tophat，因此作者推薦採用Hisat進行序列比對。Hisat對參考基因組建立索引的命令為：$HISAT_HOME/hisat-build reference_genome db_name利用Hisat進行比對的命令為：$HISAT_HOME/hisat –p 4 –x db_name -1 left.fq -2 right.fq –S output.sam如果構建的是單端reads的庫，則比對命令為：$HISAT_HOME/hisat –p 4 –x db_name –U single.fq –S output.sam輸出的比對格式為SAM格式，轉換為BAM格式命令為：$SAMTOOLS_HOME/samtools view -@ 8 –bS output.sam > output.bam得到BAM文件格式的比對結果後，就可以進行轉錄本的表達量分析。評價基因的表達量時，一種是計算轉錄本的FPKM（fragments per kilobase of transcript per million fragments）值；另一種是計算轉錄本比對上的reads的數目，通常情況下選取比對到最好的且唯一位置的reads。兩種方式的代表軟體分別有Cufflinks（Trapnell et al.，2012）、stringTie（Pertea et al.，2015）和Htseq（Anders et al.，2015）。相對於Cufflinks，stringTie的組裝速度更快，內存消耗也更小。stringTie的命令如下：$STRINGTIE_HOME/stringtie bamfile –G guide_gff -e –o out_gtfHtseq的命令如下：$HTSEQ_HOME/htseq-count -s no -i ID -t gene -q samfile gffFile > output.txt最後根據泊松分布或負二項分布模型，判斷基因表達是否在統計學水平上有差異，根據設定的差異基因篩選標準，篩選出滿足條件的差異基因列表。這一過程可以利用Cuffdiff（Trapnell et al.，2012）、Ballgown（Frazee et al.，2015）、DEseq（Anders and Huber，2010）等軟體進行差異表達基因的統計。以Cuffdiff為例，命令如下：$CUFFLINKS_HOME/cuffdiff –p 1 –u merged.gtf –b genome.fa –L S1，S2 –o S1_S2.diff S1.bam S2.bam轉錄組統計篩選差異基因時，通常採取差異比值大於2倍，校正後顯著性P值小於等於0.001。（四）可變剪切分析RNA-seq測序不僅可以分析已知轉錄本的表達量，還可以用於發現新的轉錄本，檢測不同的可變剪切模式。進行Hisat、Star等序列比對軟體得到的BAM文件，可以通過Cufflink、stringTie等軟體進行局部reads組裝，發現新的轉錄本信息。在這一過程中，如果已經具有已知的基因注釋集合，那麼可以加到數據分析中去，用來優化基因的結構。在轉錄組組裝時，常常會涉及多種樣品的組裝，如在不同時間點、不同組織中取樣的樣品，或者用不同條件處理的樣品等，所以需要對初步的組裝結果進行合併。我們可以選擇Cufflink軟體包中的cuffmerge進行合併操作。cuffmerge的命令如下：$CUFFLINKS_HOME/cuffmerge –g Guide.gtf –s genome.fa –o merge.gtf assembl_gtf_list.txt關於轉錄本的描述文件GTF通常包含上萬條記錄，我們需要應用軟體工具對它進行解析。Asprofile（Florea et al.，2013）可以對不同時間點、不同組織或不同處理條件的樣品的RNA-seq數據進行可變剪切事件的提取和比較，並做轉錄本的定量分析。對於組裝的轉錄本可以通過IGV（Robinson et al.，2011）等基因組瀏覽器進行圖形化展示。（五）SNP分析轉錄組數據可以用來開發分子標記。SNP（single nucleotide ploymorphism）即單核苷酸多態性，由於其覆蓋率極高足以用於區分兩個生物樣本和定位相應的基因，是非常熱門的一類標記，經常被用來描述基因組上DNA的變異。在轉錄組測序中，由於測序區域大多位於基因區域，轉錄組的SNP分析更多的是用來區分測序樣品在功能基因區域的差異，這些位於基因區域的差異更有可能與人們感興趣的一些表型性狀相關。利用GATK（McKenna et al.，2010b）進行轉錄組的SNP分析流程如下。通過Hisat與參考基因組進行比對：$HISAT_HOME/hisat –p 4 –x db_name -1 left.fq -2 right.fq –S output.sam添加reads測序相關信息，排序，去除重複reads：java -jar AddOrReplaceReadGroups I=star_output.sam O=rg_added_sorted.bam SO= coordinate RGID=id RGLB=library RGPL=platform RGPU=machine RGSM= samplejava -jar MarkDuplicates I=rg_added_sorted.bam O=dedupped.bam CREATE_ INDEX=true VALIDATION_STRINGENCY=SILENT M=output.metrics截除intron區域reads片段如圖4.18，可以看出經過截除後，reads比對結果準確性得到提高：java-jar GenomeAnalysisTK.jar -T SplitNCigarReads -R ref.fasta -I dedupped. bam -o split.bam -rf ReassignOneMappingQuality -RMQF 255 -RMQT 60 -U ALLOW_N_CIGAR_READS這一步可以提高SNP的準確率，去除假陽性SNP。SNP的檢索及過濾命令如下。java -jar GenomeAnalysisTK.jar -T HaplotypeCaller -R ref.fasta -I input.bam -dontUseSoftClippedBases -stand_call_conf 20.0 -stand_emit_conf 20.0 -o output.vcfjava -jar GenomeAnalysisTK.jar -T VariantFiltration -R ref.fasta -V input.vcf -window 35 -cluster 3 -filterName FS -filter "FS > 30.0" -filterName QD -filter "QD < 2.0" -o output.vcf

圖4.18 GATK截除intron區域reads的前後比對結果二、轉錄組數據的功能注釋（一）InterPro簡介為了能夠進一步了解組裝後的轉錄本序列及其可能行使的功能，我們需要對這些組裝序列進行功能注釋。進行功能注釋的方法有很多種，基本上都是利用了「序列相近、功能相近」的原則，即利用序列比對軟體進行同源比對。進行功能注釋的方法有很多，可以直接利用現有的一些功能注釋較好的資料庫，如用Uniprot人工校正過的Swissprot資料庫或者模式生物已注釋的資料庫，直接進行blast比對，設定閾值，從而得到相近的功能信息。另外，也可以使用軟體及資料庫進行功能注釋，如PFAM、SAMRT、SUPERFAMILY、InterPro（Mitchell et al.，2015）等。InterPro（ebi.ac.uk/interpro/）是一個綜合性資料庫，主要用於蛋白質及基因組的分類與自動注釋。InterPro將序列按照超家族、家族、子家族等不同水平，將蛋白質序列進行分類，預測序列的功能保守域、重複序列、關鍵位點等。其相應的軟體為InterProScan（Jones et al.，2014）。下載InterProScan後需要下載相應的各種資料庫及相應軟體進行配置。InterProScan不僅可以進行序列的功能注釋，還可以進行GO、Pathway等相應的高級功能分析，因此它是一個綜合功能比較強大的、可用於批量化注釋的軟體包。（二）COG分析Orthologs是指來自於不同物種的由垂直家系（物種形成）進化而來的蛋白質，並且典型地保留著與原始蛋白質相同的功能。Paralogs是那些在一定物種中來源於基因複製的蛋白質，可能會進化出與原來功能有關的新功能。「COG」是cluster of orthologous groups of proteins（直系同源蛋白相鄰類的聚簇）的縮寫（Tatusov et al.，1997）。由於Orthologs一般保留相同的功能，因此，利用NCBI的COG資料庫，基於同源序列比對的信息，可以推斷未知序列的功能，以及是否參與特定的代謝途徑。COG用不同的大寫字母進行功能及代謝途徑的分類。圖4.19是基於COG分析後得到的功能分類圖。

圖4.19 COG功能分類樣圖（三）GO注釋「GO」是gene ontology（基因本體論）的縮寫。它給出了3種不同的定義，包括分子功能（molecular function）、生物學途徑（biological process）、細胞組分（cellular component）。這3種不同特徵的定義可以對基因產物進行全方位的描述。對於一些已知物種，如果GO官網（geneontology.org/）已經完成了該物種的蛋白質分類，可以直接檢索使用。如果從該網站檢索不到目的蛋白，可以通過序列比對，找到含有功能注釋的同源蛋白。通過Blast2GO（Conesa et al.，2005）或者InterProScan軟體可以實現這一過程。Blast2GO是一個圖形化界面的軟體，通過分步操作，可以實現從序列比對到最終功能注釋等許多分析，同時，還可以對注釋結果進行分類匯總。由於Blast2GO簡便易學，因此推薦讀者使用Blast2GO軟體進行功能分析。圖4.20為該軟體序列比對的界面圖（blast2go. com/blast2go-pro），可以選擇不同資料庫進行比對。

圖4.20 Blast2GO軟體界面（四）Pathway注釋目前與Pathway相關的資料庫有KEGG、MetaCyc、Reactome、UniPathway等。與GO注釋相似，Pathway的注釋也需要通過序列比對進行預測。能夠完成這一過程的軟體有KOBAS2.0（Xie et al.，2011）、InterProScan，Blast2GO也有相應的模塊進行Pathway的分析。KOBAS2.0既有在線的伺服器，也可以下載後構建本地資料庫進行分析。KOBAS2.0是使用Python編寫的軟體，分為兩部分。第一部分是annotate過程，即將序列與KEGG等資料庫中的蛋白質序列進行比對，從而得到對應關係，這一過程可以幫助我們得到Unigene集合相對應的Pathway的信息。第二部分是identify過程，即功能富集分析，在本節「二、（五）」中介紹。KOBAS2.0命令行如下。python annotate.py -i query.tab -s db -t blastout:tab -o annotate.out -e 1e-20 -n 1KEGG有相應的API程序（kegg.jp/kegg/rest/kegga），可以對KOBAS2.0的結果進行豐富註解，如圖4.21所示，通過其API的介面可以清楚地顯示出不同類型的蛋白質信息。

圖4.21 KEGG中氧化磷酸化途徑展示圖（五）功能富集分析GO注釋及Pathway注釋均可以進行功能富集分析。功能富集分析採用的是一種概率演算法，稱為超幾何分布（hypergeometric distribution）。通過Fisher精確檢驗（Fisher exact test），用於描述兩個集合中是否符合同一分布的概率。在轉錄組數據分析中，參考數據的集合，即背景數據，通常以所有的表達基因作為背景。而某些特定的基因，則被認為是前景數據，通常組織或者處理間所鑒定的差異基因作為前景。通過計算前景數據與背景數據在某個GO或者Pathway分類中的超幾何分布關係，可以返回該前景數據在這個分類上的顯著差異值，即P值，再經過多重校正（multiple testing），可以得到校正後的P值。一般功能富集的篩選標準為校正後P值小於0.05。校正後的P值越小，前景數據與背景數據的差異就越顯著，表明所關注的前景基因可能與該功能分類密切相關。Blast2GO及KOBAS2.0均可以進行功能富集分析。這個過程也是很簡單的。在KOBAS2.0中命令行具體如下。python identify.py -f foreground.out -b background.out -o annotate.identify三、小RNA數據分析小RNA（miRNA）來源於一段有摺疊的發卡結構的單鏈RNA序列，主要參與調節內源基因的轉錄和翻譯（Carthew and Sontheimer，2009）。（一）數據預處理因為小RNA序列長度是18～30 nt，小於高通量測序得到的reads長度，所以reads會有一段3′接頭序列。因此在處理小RNA序列時，除了常規的去污染、去低質量序列，還需要去掉接頭序列。這一過程與DNA質控軟體相同，可以採用Trimmomatic等工具進行處理。（二）序列比對miRNA測序數據有兩個特徵，一是序列讀長短，二是相同的序列重複率高。因此在做序列比對時，可以先把序列去冗餘，統計出每條序列在一個樣品中有多少次重複，再比對到參考基因組。比對時要保證完全匹配，並且如果reads可以比對到多個位置，則要保留每個位置的比對記錄。能夠完全比對到參考基因組的序列進行下一步的分析。採用的比對軟體為BWA等。（三）序列分類將比對到基因組上的reads與已知類型序列比對，可以對這些reads進行分類。與已有的基因組注釋比對，可以區分編碼和非編碼RNA；與miRBase（Griffiths-Jones et al.，2006）資料庫比對，可以得到已知miRNA；與Rfam（Nawrocki et al.，2015）資料庫及GenBank資料庫比對，可得到rRNA、tRNA、snRNA等降解片段；與Repbase（Bao et al.，2015）資料庫比對，可以找出重複序列和轉座子相關序列。在剩下的序列中篩選出18～30 nt的序列進行進一步的miRNA預測。（四）新miRNA預測通過序列分類，篩選出可能的新miRNA序列。miRNA前體的標誌性髮夾結構能夠用來預測新的miRNA。根據序列比對結果，截取附近區域的一段序列，作為miRNA的前體序列，利用miRNA預測軟體，可從計算角度判斷該序列是否是miRNA，進一步的判斷還需要實驗驗證。相應的軟體有miRPlant（An et al.，2014）、mirtools（Wu et al.，2013）等。miRPlant是一個界面友好，且不依賴第三方軟體的預測植物miRNA的軟體（圖4.22）。從拿到原始測序reads到預測出miRNA，都可通過miRPlant獨立完成。miRPlant軟體包含以下幾個步驟。

圖4.22 miRPlant軟體界面（引自miRPlant: an integrated tool for identification of plant miRNA from RNA sequencing data. BMC Bioinformatics. 2014, 15:275）1）過濾掉長度在10～23 bp之外的reads，或質量值低於設定閾值的reads。2）序列一樣的reads合併為一條reads。3）以完全匹配的方式把reads比對到參考基因組。4）通過RNA二級結構演算法確定RNA-seq reads在參考基因組的臨近區域能否形成發卡結構。5）通過miRDeep model給預測的miRNA打分，分數越高，預測的miRNA的準確性越高。軟體需要輸入adapter序列和測序reads或者bam文件，其他參數可以根據情況自行調整。adapter序列：由於miRNA本身序列比測序reads短，在測序過程中會測到接頭序列，miRPlant提供了去接頭序列的功能。adapter序列可根據實際情況選擇對應的接頭序列。軟體可以接受fq、fa和bam文件，如果提供的是fq和fa，會首先進行序列比對，再執行miRNA識別程序，如果是bam文件，會直接執行miRNA識別程序。（五）miRNA的靶基因預測通過已有的資料庫，可以獲得已知miRNA的靶基因，新的miRNA需要通過靶基因預測軟體進行預測。通過miRNA和靶基因的對應關係，我們可以對基因的調控關係進行研究。這一過程中的預測軟體有psRobot（Wu et al.，2012）、miRBase、miRNAMap（Hsu et al.，2008）和TargetScanS（Lewis et al.，2005）等。psRobot提供在線和本地兩種版本（參考omicslab.genetics.ac.cn help_stem_loop.php），功能如下。1）把小RNA比對到參考基因組，確定其在基因組上的位置。2）分析小RNA是否為重複序列及重複序列的類型。3）小RNA在小RNA合成相關蛋白複合物中的表達情況。4）小RNA前體序列的二級結構。5）預測小RNA的目標基因。6）小RNA和目標基因在擬南芥、水稻、高粱、玉米等植物中的保守性。7）降解組數據中目標基因的剪切信息。預測小RNA的目標基因：需要提供的數據包括小RNA的序列和mRNA的序列，這兩類數據都可以用已有的序列，或選擇自己提供序列。然後設定預測目標基因的相關參數，即可進行小RNA的目標基因預測。（六）表達模式分析通過統計樣品中miRNA的數目，可以獲得count-base類型的表達量。可以對樣品進行聚類分析、差異表達分析、通路及功能的富集分析，具體分析參考轉錄組分析流程。（曹英豪 李 妍）第四節 建造中等高性能計算機群系統一、高性能計算技術簡介高性能計算（HPC）是一個計算機集群系統，它通過各種互聯技術將多個計算機系統連接在一起，利用被連接的系統的綜合計算能力來處理大型數據與計算問題，通常又被稱為高性能計算集群。高性能計算方法的基本原理就是將問題分為若干部分，集群中的每台計算機（稱為節點）均可同時參與問題的解決，從而顯著縮短了計算時間。解決大型計算問題需要功能強大的計算機系統。隨著高性能計算系統的出現，這一類應用從昂貴的大型外部計算機系統演變為採用商用伺服器產品和軟體的高性能計算機集群。因此，高性能計算系統已經成為解決大型問題的計算機系統的發展方向。那麼，什麼樣的大型問題最適合使用高性能計算系統呢？一般來說，高性能計算系統是為了滿足以下需要的計算系統：①能夠突破性能極限的計算；②單個高端計算機系統不能滿足其需求的計算；③需要通過專門的程序優化最大限度提高系統的I/O、計算和數據傳送性能的計算。二、如何有效地構建高性能計算集群一套成熟的、高效的高性能計算集群，應由以下幾部分構成：節點、存儲系統、作業調度系統、高效計算網路、加速卡等。下面我們對這五項分別進行闡述。（一）節點高性能計算是由多個節點組成的，每個節點是一個獨立的管理或計算單元，能夠獨立參與任務的計算，節點數量及節點性能是反映一個高性能計算集群計算能力的基本指標。（二）存儲系統高性能計算中的存儲系統有其獨特性。由於高性能計算集群具有節點眾多、需要各節點協同完成同一任務等特性，因此要求接入高性能計算的存儲系統應具有高並發、高I/O、高處理效率、統一命名空間等特性。（三）作業調度系統由眾多計算、管理、存儲等節點共同參與的高性能計算集群系統，必須有一套高效快捷的系統把這些資源管理起來，它就是作業調度系統，這套系統是高性能計算集群的重中之重。作業調度系統的主要功能是根據作業控制塊中的信息，審查系統能否滿足用戶作業的資源需求，以及按照一定的演算法，從外存的後備隊列中選取某些作業調入內存，並為它們創建進程、分配必要的資源。然後再將新創建的進程插入就緒隊列，準備執行。因此，有時也把作業調度稱為接納調度。（四）高效計算網路高性能計算集群中一大類任務是mpi類計算任務，該任務需要多個節點共同協作來完成。一個複雜的作業由作業調度系統拆解成若干個計算單元，再將這些計算單元分配到不同節點的不同計算核心中進行並行計算。由於並行計算是具有結果依賴性的，因此在每一步或每幾步計算之後，所有節點的計算核心要進行數據交互。這就需要我們的計算網路具有高帶寬、低延遲特性。假設我們的並行顆粒度相對較大，儘管單位任務被拆解得非常小，但由於一次計算可能要經曆數百萬次的數據交互，每次數據交互如果有0.1 s的延遲，那麼數百萬次的數據交互所產生的等待時間是讓人無法忍受的，甚至會超過任務拆解前所需要的計算時間。所以高效的計算網路，是保證並行計算有效運行的基本前提。（五）加速卡高性能計算集群中的主要計算資源是CPU，但一些特殊的場景或一些特殊的演算法可以通過GPU或MIC來進行加速，從而達到提高計算效率的目的。三、建造中等測序中心所需高性能計算集群的硬體配置隨著測序技術的飛速發展，單位時間內的數據產出量呈幾何倍數增長，從而推進了測序中心中高性能計算集群技術的發展。下面我們根據目前中等測序中心的數據產出及處理能力，給出一個相對的高性能計算集群的基本配置需求。（一）系統架構整套集群採用高性能計算系統架構。配有作業調度系統，各節點採用Linux CentOS6操作系統，兩個管理節點配置High Available，4個登錄節點配置Load Balance，從而達到整個集群系統高度可用且可負載均衡。（二）計算能力總計算能力不低於15 Tflops/s，其中應按照1︰4的計算核心比例來配置胖節點與普通計算節點。胖節點的總內存數應不低於1 TB，CPU核心數應大於80核。計算節點總內存數應不低於196 GB，CPU核心數應大於24核。（三）存儲需求應採用並行存儲系統，I/O節點數目不少於10個，MDS節點應配置HA。裸容量不低於1PB，總帶寬不小於4 GB/S，可負載不小於3 W的IOPS。該存儲支持並行文件系統，且在集群中以統一命名空間訪問。（四）計算網路計算網路應採用56 Gb/s的無限帶寬技術（infiniband）高速計算網路，保證整套系統的高帶寬及低延時。四、高性能計算集群的工作環境要求伺服器存儲設備需要放置在專用的中心機房中使用。為了保證高性能計算集群的高效運行，高性能計算集群所依賴的中心機房尤為重要。下面簡單介紹高性能計算中心機房的建設標準及工作環境要求。機房建設是建築智能化系統的一個重要部分。機房建設涵蓋了建築裝修、供電、照明、防雷、接地、UPS不間斷電源、精密空調、環境監測、火災報警及滅火、門禁、防盜、閉路監視、綜合布線和系統集成等技術。（一）高性能計算中心機房建設思想整體機房建設：將機房設備、監控設備、強弱電系統、數據/非數據設備等作為一個完整的系統考慮，盡量發揮各子系統的聯動、互動作用。可管理、可擴展：現代機房建設已不僅僅是功能上的要求，而且要具有良好的可管理性，為用戶提供友善的管理界面，同時要保證容量、性能的可擴展性，以保護用戶投資。高質量項目管理：機房建設是一項專業化的綜合性工程。要求對裝修、配電、空調、新排風、監控、門禁、消防、防雷、防過壓、接地、綜合布線和網路等各個子系統的建設規劃、方案設計、施工安裝等過程進行嚴密的統籌管理，以保證工程的質量和周期。（二）機房建設標準機房建設應遵循如下標準。國家標準《電子計算機機房設計規範》（GB50174-93）國家標準《計算站場地技術條件》（GB2887-89）國家標準《電子計算機機房施工及驗收規範》（SJ/T30003-93）國家標準《計算機機房活動地板的技術要求》（GB6650-86）國家標準《計算站場地安全技術》（GB9361-88）國家標準《電氣裝置安裝工程接地裝置施工及驗收規範》（GB50169-92）公共安全行業標準《安全防範工程程序與要求》（GA/T75-94）《中華人民共和國計算機信息系統安全保護條例》《工業企業照明設計標準》（GB50034-92）《室內裝飾工程質量規範》（QB 1838-93）（三）機房建設方案機房建設方案的設計應根據用戶提出的技術要求，對機房建設的建築物進行實地勘查，依據國家有關標準和規範，結合所建各種系統運行特點進行總體設計。總體機房建設方案以業務完善技術規範，安全可靠為主，確保系統安全可靠地運行。在選材投資方面根據功能及設備要求區別對待，並滿足用戶的特殊要求，做到投資有重點，保證機房場地工作人員的身心健康，延長系統的使用壽命。機房建設的工作就是圍繞這個根本任務，通過採用優質產品和先進工藝把上述設計思想有機地結合起來，為機房裡的設備和工作人員創造一個安全、可靠、美觀、舒適的工作場地。一個全面的機房建設應包括以下幾個方面：①機房裝修；②供電系統；③空調系統；④門禁系統；⑤監控系統；⑥消防報警系統；⑦防雷接地系統。1.機房裝修機房建設中的機房裝修主要包括吊頂，隔斷牆，門、窗、牆壁裝修，地面、活動地板的施工驗收及其他室內作業。機房裝修是整個機房建設的基礎，是機房環境建設的重要環節。機房必須具備防塵、防潮、抗靜電、阻燃、絕緣、隔熱、降雜訊的物理環境，機房功能區域分隔要清晰明了、便於識別和維護。機房裝修作業應符合《裝飾工程施工及驗收規範》、《地面與樓面工程施工及驗收規範》、《木結構工程施工及驗收規範》及《鋼結構工程施工質量驗收規範》的有關規定。在機房建設過程中應保證現場、材料和設備的清潔。隱蔽工程（如地板下、吊頂上、假牆、夾層內）在封口前必須先除塵、清潔處理，暗處表層應能保持長期不起塵、不起皮和不龜裂。機房所有管線穿牆處的裁口必須做防塵處理，必須對縫隙採用密封材料填堵。在裱糊、粘接貼面及進行其他塗復施工時，其環境條件應符合材料說明書的規定。機房裝修材料應盡量選擇無毒、無刺激性、難燃、阻燃的材料，否則應儘可能塗防火塗料。2.供電系統1）保證機房核心設備（伺服器、網路設備及輔助設備）安全穩定運行。核心設備供電系統必須達到一類供電標準，即必須建立不間斷供電系統。市電停電時，不間斷電源系統主機電源實際輸出功率宜大於後端負載1.5倍，滿負荷運轉時間不得少於120 min。2）信息系統設備供電系統必須與動力、照明系統分開。供電系統要求：頻率，50 Hz；電壓，380 V/220 V；相數，三相五線或三相四線制/單相三線制；穩態電壓偏移範圍220 V（–10%～+10%）；穩態頻率偏移範圍，50Hz（–0.5%～+0.5%）。3）電力布線要求：機房UPS電源要採用獨立雙迴路供電，輸入電流應符合UPS輸入端電流要求；靜電地板下的供電線路置於管內，分支到各用電區域，向各個用電插座分配電力，防止外界電磁干擾系統設備；線路上要有標籤說明去向及功能。3.空調系統1）通過機房空調系統保持機房內相對穩定的溫度和濕度，使機房內的各類設備保持良好的運行環境，確保系統可靠、穩定運行。機房系統要求：全年溫度，18～25℃；相對濕度，35%～65%；溫度變化率，<10℃/h；並不得結露。2）機房要有防水害措施，確保機房安全運行（機房一般都配備恆溫濕裝置，所以在一般情況下禁止使用冷、暖氣系統，如已使用則必須對系統給排水管道採取嚴格的防漏及補救措施）。3）根據《計算站場地技術條件》（GB2887-89）機房技術要求，按A級設計，溫度T=（23±2）℃，相對濕度=55%±5%，夏季取上限，冬季取下限。氣流組織採用下送風、上迴風，即抗靜電活動地板靜壓箱送風，吊頂天花微孔板迴風。新風量設計取總風量的10%，中低度過濾，新風與迴風混合後，進入空調設備處理，提高控制精度，節省投資，方便管理。4.門禁系統機房建設中的門禁管理系統的主要目的是保證重要區域設備和資料的安全，便於人員的合理流動，對進入這些重要區域的人員實行各種方式的門禁管理，以便限制人員隨意進出；卡片最好採用現在流行的感應式卡片；卡出入系統首先應具有許可權設置的功能，即每張卡可進出的時間、可進出哪道門，不同的卡片持有者應有不同的許可權；每次有效的進入都應存檔或統計；應有完善的密碼系統，即對系統的更改，不同的操作者應有不同的許可權；電鎖應採用安全可靠的產品，有電閉鎖或無電閉鎖根據用戶要求可調；緊急情況下或電鎖出現故障的情況下應有應急鑰匙可將門打開；門禁系統最好採用計算機控制系統；全套系統最好有備用電源。5.監控系統機房中有大量的伺服器及機櫃、機架。由於這些機櫃及機架一般比較高，監控的死角比較多，因此在電視監控布點時主要考慮各個出入口，每一排機櫃之間安裝攝像機。如果在各出入口的空間比較大，可考慮採用帶變焦的攝像機，在每一排的機櫃之間，根據監視距離，配定焦攝像機即可。如果機房有多個房間的話，可考慮在UPS房和控制機房內安裝攝像機。機房監控圖像信號應保持24 h錄像，錄像方式可採用硬碟錄像，也可採用傳統的錄像系統。閉路電視控制系統最好有視頻動態報警功能。同時如果具有視頻遠程傳輸功能，即通過互聯網、ISDN、區域網或電話線將監視信號傳輸到遠程客戶指定的地方，在使用時將會更加方便。6.消防報警系統機房內有許多重要設備，其價值較高，分布密集，且需要24 h不間斷運行，對消防要求較高，一旦機房出現火災，為保證設備安全，需要消防系統在第一時間內自動啟動進行滅火，通知遠程管理人員減少損失。機房的物理環境：機房的結構、材料、配置設施必須滿足保溫、隔熱、防火等要求。機房應有溫感、煙感、報警器等裝置和消防設備設施，必須採用氣體滅火劑。7.防雷接地系統為防止機房設備的損壞和數據的丟失，機房建設中機房的防雷接地尤其重要。按國家建築物防雷設計規範，本機房建設方案對機房電氣電子設備的外殼、金屬件等實行等電位連接，並在低壓配電電源電纜進線輸入端加裝電源防雷器。機房防雷接地系統注意以下兩點。1）機房的接地系統必須安裝室外的獨立接地體；直流地、防靜電地採用獨立接地；交流工作地、安全保護地採用電力系統接地；不得共用接地線纜，所有機櫃必須接地。2）機房防雷系統的外部防護主要由建築物自身防雷系統來承擔；由室外直接接入機房金屬信息線纜，必須作防浪涌處理；所有弱電線纜不裸露於外部環境；弱電橋架使用扁銅軟線帶跨接，進行可靠接地；機房電源系統至少進行二極防浪涌處理；重要負載末端防浪涌處理。機房建設除了上述系統外，還有穿插於整個機房建設的綜合布線系統，機房綜合布線系統應滿足主幹線路有冗餘，機房布線系統中所有的電纜、光纜、信息模塊、接插件、配線架、機櫃等在其被安裝的場所均要容易被識別，線纜布設整齊，布線中的每根電纜、光纜、信息模塊、配線架和端點要指定統一的標誌符，電纜在兩端要有標註，保證維修方便、操作靈活。五、下機數據的自動化傳輸以HiSeq2500為例，測序儀所配置主機的存儲空間及存儲效率有限，且最終數據要放到高性能計算集群中進行分析。所以我們有必要對下機數據進行自動化傳輸配置，保證測序下機數據自動傳輸到高性能計算集群中，方便後續的數據處理。針對類似於HiSeq系列的測序儀，控制器操作系統為Windows系統，我們可以採用CIFS或ISCSI的方式將存儲系統掛載到測序儀上。將測序數據的生成目錄更改為此動態存儲，即可完成下機數據的自動傳輸。針對類似於Ion Torrent系列的測序儀，控制器操作系統為類Linux系統，我們可以採用NFS或ISCSI的方式將存儲系統掛載到測序儀上。將測序數據的生成目錄更改為此動態存儲，即可完成下機數據的自動傳輸。在數據處理段，應將該共享卷設置為只讀。下機數據的預處理不需要對原始下機數據進行任何寫操作，從而保證數據的讀一致性。（楊 鑫）參 考 文 獻Abe A, Kosugi S, Yoshida K, et al. 2012. Genome sequencing reveals agronomically important loci in rice using MutMap. Nature Biotechnology, 30(2): 174-178.An J, Lai J, Sajjanhar A, et al. 2014. miRPlant: an integrated tool for identification of plant miRNA from RNA sequencing data. BMC Bioinformatics, 15: 275.Anders S, Huber W. 2010. Differential expression analysis for sequence count data. Genome Biol, 11(10): R106.Anders S, Pyl P T, Huber W. 2015. HTSeq--a Python framework to work with high-throughput sequencing data. Bioinformatics, 31(2): 166-169.Bao W, Kojima K K, Kohany O. 2015. Repbase Update, a database of repetitive elements in eukaryotic genomes. Mob DNA, 6: 11.Belfield E J, Gan X, Mithani A, et al. 2012. Genome-wide analysis of mutations in mutant lineages selected following fast-neutron irradiation mutagenesis ofArabidopsis thaliana. Genome Research, 22(7): 1306-1315.Bolger A M, Lohse M, Usadel B. 2014. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics: btu170.Butler J, MacCallum I, Kleber M, et al. 2008. ALLPATHS: de novo assembly of whole-genome shotgun microreads. Genome Res, 18(5): 810-820.Carthew R W, Sontheimer E J. 2009. Origins and Mechanisms of miRNAs and siRNAs. Cell, 136(4): 642-655.Chen W, Gao Y, Xie W, et al. 2014a. Genome-wide association analyses provide genetic and biochemical insights into natural variation in rice metabolism. Nature Genetics, 46(7): 714-721.Chen Z, Wang B, Dong X, et al. 2014b. An ultra-high density bin-map for rapid QTL mapping for tassel and ear architecture in a large F2 maize population. BMC Genomics, 15(1): 433.Chikhi R, Medvedev P. 2013. Informed and automated k-mer size selection for genome assembly. Bioinformatics: btt310.Conesa A, Gotz S, Garcia-Gomez J M, et al. 2005. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics, 21(18): 3674-3676.Danecek P, Auton A, Abecasis G, et al. 2011. The variant call format and VCFtools. Bioinformatics, 27(15): 2156-2158.Dobin A, Davis C A, Schlesinger F, et al. 2013. STAR: ultrafast universal RNA-seq aligner. Bioinformatics, 29(1): 15-21.Florea L, Song L, Salzberg S L. 2013. Thousands of exon skipping events differentiate among splicing patterns in sixteen human tissues. F1000Res, 2: 188.Frazee A C, Pertea G, Jaffe A E, et al. 2015. Ballgown bridges the gap between transcriptome assembly and expression analysis. Nat Biotechnol, 33(3): 243-246.Gordon A, Hannon G. 2010. Fastx-toolkit. FASTQ/A short-reads preprocessing tools. unpublished)http: //hannonlab. cshl. edu/fastx_toolkit [2015-8-12].Griffiths-Jones S, Grocock R J, van Dongen S, et al. 2006. miRBase: microRNA sequences, targets and gene nomenclature. Nucleic Acids Res, 34(Database issue): D140-144.Haas B J, Papanicolaou A, Yassour M, et al. 2013. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc, 8(8): 1494-1512.Hsu S D, Chu C H, Tsou A P, et al. 2008. miRNAMap 2.0: genomic maps of microRNAs in metazoan genomes. Nucleic Acids Res, 36(Database issue): D165-169.Huang X, Feng Q, Qian Q, et al. 2009. High-throughput genotyping by whole-genome resequencing. Genome Research, 19(6): 1068-1076.Jia J, Zhao S, Kong X, et al. 2013. Aegilops tauschii draft genome sequence reveals a gene repertoire for wheat adaptation. Nature, 496(7443): 91-95.Jiang Y H, Yuen R K, Jin X, et al. 2013. Detection of clinically relevant genetic variants in autism spectrum disorder by whole-genome sequencing. The American Journal of Human Genetics, 93(2): 249-263.Jones P, Binns D, Chang H Y, et al. 2014. InterProScan 5: genome-scale protein function classification. Bioinformatics, 30(9): 1236-1240.Kim D, Langmead B, Salzberg S L. 2015. HISAT: a fast spliced aligner with low memory requirements. Nat Methods, 12(4): 357-360.Klein R J, Zeiss C, Chew E Y, et al. 2005. Complement factor H polymorphism in age-related macular degeneration. Science, 308(5720): 385-389.Lewis B P, Burge C B, Bartel D P. 2005. Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets. Cell, 120(1): 15-20.Li H, Durbin R. 2009. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics, 25(14): 1754-1760.Li H, Handsaker B, Wysoker A, et al. 2009. The sequence alignment/map format and SAMtools. Bioinformatics, 25(16): 2078-2079.Li R, Fan W, Tian G, et al. 2010. The sequence and de novo assembly of the giant panda genome. Nature, 463(7279): 311-317.Li R, Li Y, Kristiansen K, et al. 2008. SOAP: short oligonucleotide alignment program. Bioinformatics, 24(5): 713-714.Ling H Q, Zhao S, Liu D, et al. 2013. Draft genome of the wheat A-genome progenitorTriticum urartu. Nature, 496(7443): 87-90.Lynam C P, Gibbons M J, Axelsen B E, et al. 2006. Jellyfish overtake fish in a heavily fished ecosystem. Current Biology, 16(13): R492-R493.McKenna A, Hanna M, Banks E, et al. 2010. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome research, 20(9): 1297-1303.Mitchell A, Chang H Y, Daugherty L, et al. 2015. The InterPro protein families database: the classification resource after 15 years. Nucleic Acids Res, 43(Database issue): D213-221.Morris G P, Ramu P, Deshpande S P, et al. 2013. Population genomic and genome-wide association studies of agroclimatic traits in sorghum. Proceedings of the National Academy of Sciences, 110(2): 453-458.Nawrocki E P, Burge S W, Bateman A, et al. 2015. Rfam 12.0: updates to the RNA families database. Nucleic Acids Res, 43(Database issue): D130-137.Pertea M, Pertea G M, Antonescu C M, et al. 2015. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol, 33(3): 290-295.Purcell S, Neale B K, Todd-Brown, et al. 2007. PLINK: a tool set for whole-genome association and population-based linkage analyses. The American Journal of Human Genetics, 81(3): 559-575.Robertson G, Schein J, Chiu R, et al. 2010. De novo assembly and analysis of RNA-seq data. Nat Methods, 7(11): 909-912.Robinson J T, Thorvaldsdottir H, Winckler W, et al. 2011. Integrative genomics viewer. Nat Biotechnol, 29(1): 24-26.Simpson J T, Wong K, Jackman S D, et al. 2009. ABySS: a parallel assembler for short read sequence data. Genome Res, 19(6): 1117-1123.Surget-Groba Y, Montoya-Burgos J I. 2010. Optimization of de novo transcriptome assembly from next-generation sequencing data. Genome Res, 20(10): 1432-1440.Takagi H, Abe A, Yoshida K, et al. 2013. QTL-seq: rapid mapping of quantitative trait loci in rice by whole genome resequencing of DNA from two bulked populations. The Plant Journal, 74(1): 174-183.Tatusov R L, Koonin E V, Lipman D J. 1997. A genomic perspective on protein families. Science, 278(5338): 631-637.Trapnell C, Pachter L, Salzberg S L. 2009. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics, 25(9): 1105-1111.Trapnell C, Roberts A, Goff L, et al. 2012. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc, 7(3): 562-578.Wu H J, Ma Y K, Chen T, et al. 2012. PsRobot: a web-based plant small RNA meta-analysis toolbox. Nucleic Acids Res, 40(Web Server issue): W22-28.Wu J, Liu Q, Wang X, et al. 2013. mirTools 2.0 for non-coding RNA discovery, profiling, and functional annotation based on high-throughput sequencing. RNA Biol, 10(7): 1087-1092.Xie C, Mao X, Huang J, et al. 2011. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res, 39(Web Server issue): W316-322.Xie Y, Wu G, Tang J, et al. 2014. SOAPdenovo-Trans: de novo transcriptome assembly with short RNA-Seq reads. Bioinformatics, 30(12): 1660-1666.Yoshida S, Ishida J K, Kamal N M, et al. 2010. A full-length enriched cDNA library and expressed sequence tag analysis of the parasitic weed,Striga hermonthica. BMC plant biology, 10(1): 55.Zerbino D R. 2010. Using the Velvet de novo assembler for short-read sequencing technologies. Curr Protoc Bioinformatics Chapter, 11: Unit 11, 15.Zhang G, Fang X, Guo X, et al. 2012. The oyster genome reveals stress adaptation and complexity of shell formation. Nature, 490(7418): 49-54.Zhou X, Stephens M. 2012. Genome-wide efficient mixed-model analysis for association studies. Nature Genetics, 44(7): 821-824.Zhou Z, Jiang Y, Wang Z, et al. 2015. Resequencing 302 wild and cultivated accessions identifies genes related to domestication and improvement in soybean. Nature Biotechnology, 33(4): 408-414.Zimin A V, Mar?ais G, Puiu D, et al. 2013. The MaSuRCA genome assembler. Bioinformatics, 29(21): 2669-2677.
推薦閱讀：
相关文章