如何對基因組序列進行註釋

基因組組裝完成後，或者是完成了草圖，就不可避免遇到一個問題，需要對基因組序列進行註釋。註釋之前首先得構建基因模型，有三種策略：

從頭註釋(de novo prediction)：通過已有的概率模型來預測基因結構，在預測剪下位點和UTR區準確性較低
同源預測(homology-based prediction)：有一些基因蛋白在相近物種間的保守型搞，所以可以使用已有的高質量近緣物種註釋資訊通過序列聯配的方式確定外顯子邊界和剪下位點
基於轉錄組預測(transcriptome-based prediction)：通過物種的RNA-seq資料輔助註釋，能夠較為準確的確定剪下位點和外顯子區域。

每一種方法都有自己的優缺點，所以最後需要用EvidenceModeler(EVM)和GLEAN工具進行整合，合併成完整的基因結構。基於可靠的基因結構，後續可才是功能註釋，蛋白功能域註釋，基因本體論註釋，通路註釋等。

那麼基因註釋重要嗎？可以說是非常重要了，尤其是高通量測序非常便宜的現在。你可以花不到一萬的價格對600M的物種進行100X的普通文庫測序，然後拼接出草圖。但是這個草圖的價值還需要你進行註釋後才能顯現出來。有可能你和諾貝爾獎就差一個註釋的基因組。

基因註釋的重要性

從案例中學習套路

陸地棉基因組註釋

文章標題為“Sequencing of allotetraploid cotton (Gossypium hirsutum L. acc. TM-1) provides a resource for fiber improvement”.

同源註釋：從Phytozome上下載了7個植物的基因組蛋白序列(Arabidopsis thaliana, Carica papaya, Glycine max, G. raimondii, Populus trichocarpa, Theobroma cacao and Vitis vinifera), 使用 TblastN

將蛋白序列比對到組裝序列上，E-value的閾值為1e-5. 將不同蛋白的BLAST的hits用 Solar 軟體進行合併。GeneWise 根據每個BLAST hit的對應基因區域預測完整的基因結構。

從頭預測：先得構建repeat-mask genome，在這個基礎上就用 August, Genescan, GlimmerHMM, Geneid 和 SNAP 預測編碼區

轉錄組預測：用Tophat將RNA-seq資料比對到組裝序列上，然後用cufflinks組裝轉錄本形成基因模型。

綜上，使用 EvidenceModeler(EVM) 將上面的結果組裝成非冗餘的基因結構。進一步根據Cscore > 0.5，peptide coverage > 0.5 和CDS overlaping with TE進行篩選。還有過濾掉超過30%編碼區被Pfam或Interprot TE domain的註釋的基因模型。

這些基因模型使用BLASTP進行功能註釋，所用資料庫為SWiss-Prot和TrEMBL.蛋白功能使用InterProScan和HMMER註釋，資料庫為InterPro和Pfam。GO註釋則是直接僱傭InterPro和Pfam註釋得到的對應entry。通路註釋使用KEGG資料庫。

Cardamine hirsuta基因組註釋

文章標題為“The Cardamine hirsuta genome offers insight into the evolution of morphological diversity”。

同源註釋：使用 GenomeThreader 以擬南芥為剪下模型，以及PlantsGDB resourc上 Brassica rapa (v1.1), A. thaliana(TAIR10), A. lyrata (v6), tomato (v3.6), poplar (v2) 和 A. thaliana (version PUT-169), B. napus (version PUT-172) EST assemblies 的完整的代表性蛋白集。

轉錄本預測：將 C. hirsuta RNA-seq資料比對到基因序列，然後用cufflinks拼接

從頭預測：轉錄本預測得到的潛在蛋白編碼轉錄本使用網頁工具 ORFpredictor 進行預測，同時用 blastx 和 A. thalina 進行比較，選擇90%序列相似度和最高5%長度差異的部分從而保證保留完整的編碼框(有啟動子和終止子)。這些基因模型根據相互之間的相似度和重疊度進行聚類，高度相似(>95)從聚類中剔除，保證非冗餘訓練集。為了訓練gene finder, 它們選隨機選取了2000個位點，20%是單個外顯子基因。從頭預測工具為 August , GlimmerHMM, Geneid 和 SNAP . 此外還用了Fgenesh+, 以雙子葉特異矩陣為引數進行預測。

最後使用JIGSAW演算法根據以上結果進行訓練，隨後再次用JIGSAW對每個基因模型計算統計學權重。

可變剪下模型則是基於苗、葉、花和果實的RNA-seq比對組裝結果。

小結

舉的2個例子都是植物，主要是植物基因組不僅是組裝，註釋都是一大難題。因為植物基因組有大量的重複區，假基因，還有很多新的蛋白編碼基因和非編碼基因，比如說玉米基因組80%以上都是重複區域。然後當我檢索這兩篇文章所用工具的時候，我不經意或者說不可避免就遇到了這個網站 http://www.plantgdb.org/ , 一個整合植物基因組學工具和資源的網站，但是這個網站似乎2年沒有更新了。當然這個網站也挺不錯,http://bioservices.usd.edu/gsap.html, 他給出了一套完整的註釋流程以及每一步的輸入和輸出情況。

此外，2017年在《Briefings in Bioinformatics》發表的”Plant genome and transcriptome annotations: from misconceptions to simple solution” 則是從五個角度對植物基因組註釋做了很完整的總結

植物科學的常見本體
功能註釋的常用資料庫和資源
已註釋的植物基因組意味著什麼
一個自動化註釋流程
一個參考流程圖，用來說明使用公用資料庫註釋植物基因組/轉錄組的常規步驟

註釋流程圖

以上，通過套路我們對整個基因組註釋有一個大概的瞭解，後續就需要通過實際操作來理解細節。

基因組註釋

當我們談到基因註釋的時候，我們通常認為註釋是指“對基因功能的描述”，比如說A基因在細胞的那個部分，通過招募B來調控C，從而引起病變。但是基因結構也是註釋的一種形式，而且是先決條件，也就是在看似隨機的ATCG的鹼基排列中找到特殊的部分，而這些特殊的區域有著不一樣的功能。

gene structure

在正式啟動基因組註釋專案之前，需要先檢查組裝是否合格，比如contig N50的長度是否大於基因的平均長度，使用BUSCO/CEGMA檢查基因的完整性，如果不滿足要求，可能輸出結果中大部分的contig中都不存在一個完整的基因結構。當組裝得到的contig符合要求時，就可以開始基因組註釋環節，這一步分為三步：基因結構預測，基因功能註釋，視覺化和質控。

基因組結構註釋

基因結構註釋應是功能註釋的先決條件，完整的真核生物基因組註釋流程需要如下步驟：

必要的基因組重複序列遮蔽
從頭尋找基因, 可用工具為: GeneMarkHMM, FGENESH, Augustus, SNAP, GlimmerHMM, Genscan
同源蛋白預測, 內含子分析: GeneWIse, Exonerate, GenomeThreader
將EST序列，全長cDNA序列和Trinity/Cufflinks/Stringtie組裝的轉錄組和基因組聯配
如果第4步用到了多個數據來源，使用PASA基於重疊情況進行聯配
使用EvidenceModler根據上述結果進行整合
使用PASA更新EVM的一致性預測，增加UTR註釋和可變剪下註釋
必要的人工檢查

基本上是套路化的分析流程，也就有一些工具通過整合幾步開發了流程管理工具，比如說BRAKER結合了GeneMark和Augustus，MAKER2整合了SNAP,Exonerate，雖然BRAKER說自己的效果比MAKER2好，但是用的人似乎不多，根據web of knowledge統計，兩者的引用率分別是44,283, 當然BRAKER是2016，MAKER2是2011，後者在時間上有優勢。

這裡準備先按部就班的按照流程進行註釋，所用的資料是 Cardamine hirsuta , 資料下載方式如下

# Cardamine hirsutat基因組資料
mkdir chi_annotation && cd chi_annotation
wget http://chi.mpipz.mpg.de/download/sequences/chi_v1.fa
cat chi_v1.fa | tr 'atcg' 'ATCG' > chi_unmasked.fa
# 註釋結果
wget http://chi.mpipz.mpg.de/download/annotations/carhr38.gff
# Cardamine hirsutat轉錄組資料
mkdir rna-seq && cd rna-seq
wget -4 -q -A '*.fastq.gz' -np -nd -r 2 http://chi.mpipz.mpg.de/download/fruit_rnaseq/cardamine_hirsuta/ &
wget -4 -q -A '*.fastq.gz' -np -nd -r 2 http://chi.mpipz.mpg.de/download/leaf_rnaseq/cardamine_hirsuta/ &

軟體安裝不在正文中出現，會放在附錄中，除了某些特別複雜的軟體。

01-重複序列遮蔽

重複遮蔽：真核生物的基因組存在大量的重複序列，植物基因組的重複序列甚至可以高達80%。儘管重複序列對維持染色體的空間結構、基因的表達調控、遺傳重組等都具有重要作用，但是卻會導致BLAST的結果出現大量假陽性，增加基因結構的預測的計算壓力甚至影響註釋正確性。基因組中的重複按照序列特徵可以分為兩類：串聯重複(tandem repeats)和散在重複(interspersed repeats).

人類中的重複序列劃分

鑑定基因組重複區域的方法有兩種：一種基於文庫(library)的同源(homology)方法，該文庫收集了其他物種的某一種重複的一致性序列，通過相似性來鑑定重複；另一種是從頭預測(de novo)，將序列和自己比較或者是高頻K-mer來鑑定重複。

目前重複序列註釋主要軟體就是RepeatMasker和RepeatModel。這裡要注意分析的fasta的ID不能過長，不然會報錯。如果序列ID過長可以使用bioawk進行轉換，後續用到RepatModel不支援多行存放序列的fasta格式。

直接使用同源註釋工具RepeatMasker尋找重複序列：

mkdir 00-RepeatMask
~/opt/biosoft/RepeatMasker/RepeatMasker -e ncbi -species arabidopsis -pa 40 -gff -dir 00-RepeatMask/ chi_unmasked.fa
# -e ncbi
# -species 選擇物種 用~/opt/biosoft/RepeatMasker/util/queryRepeatDatabase.pl -tree 瞭解
# -lib 增加額外資料庫,
# -pa 平行計算
# -gff 輸出gff註釋
# -dir 輸出路徑
# annotation with the library produced by RepeatModel

輸出結果中主要關注如下三個(其中xxx表示一類檔名)

xxx.fa.masked, 將重複序列用N代替
xxx.fa.out.gff, 以gff2形式存放重複序列出現的位置
xxx.fa.tbl, 該檔案記錄著分類資訊

cat 00-RepeatMask/chi_unmasked.fa.tbl
==================================================
file name: chi_unmasked.fa
sequences:           624
total length:  198654690 bp  (191241357 bp excl N/X-runs)
GC level:         35.24 %
bases masked:   35410625 bp ( 17.83 %)
==================================================

也就是說該物種198M中有將近18%的重複序列，作為參考，擬南芥125Mb 14%重複序列, 水稻389M，36%重複，人類基因組是3G，50%左右的重複序列。

使用最後的chi_unmasked.fa.masked用於下一步的基因結構預測。

注：當然也可以用RepeatModel進行從頭預測，得到的預測結果後續可以整合到RepeatMasker

# de novo predict
~/opt/biosoft/RepeatModeler-open-1.0.11/BuildDatabase -name test -engine ncbi output.fa
~/opt/biosoft/RepeatModeler-open-1.0.11/RepeatModeler -database test

這一步速度極其慢，由於我們的目的只是獲取遮蔽後序列降低後續從頭預測的壓力，所以可以先不做這一步。在後續分析重複序列在基因組進化上的作用時可以做這一步。下

如果從頭預測的結果與同源預測的結果有30%以上的overlap，並且分類不一致，會把從頭預測的結果過濾掉。從頭預測與同源預測結果有overlap，但是分類一致的，都會保留。但是統計的時候不會重複統計。

02-從頭(ab initio)預測基因

基於已有模型或無監督訓練

目前的從頭預測軟體大多是基於HMM(隱馬爾科夫鏈)和貝葉斯理論，通過已有物種的註釋資訊對軟體進行訓練，從訓練結果中去推斷一段基因序列中可能的結構，在這方面做的最好的工具是AUGUSTUS 它可以僅使用序列資訊進行預測，也可以整合EST, cDNA, RNA-seq資料作為先驗模型進行預測。

AUGUSTUS的無root安裝比較麻煩，我折騰了好幾天最後卒，不過辛虧有bioconda，conda create -n annotation augustus=3.3.

它的使用看起來很簡單，我們可以嘗試使用一段擬南芥已知的基因序列讓其預測，比如前8k序列

seqkit faidx TAIR10.fa Chr1:1-8000 > test.fa
augustus --speices=arabidopsis test.fa > test.gff

如果僅僅看兩者的CDS區，結果完全一致，相當於看過一遍參考答案去做題目，題目都做對了。

注:已經被訓練的物種資訊可以用augustus --species=help檢視。

結果比較

在不使用RNA-seq資料的情況下，可以基於擬南芥的訓練模型進行預測，採用下面的方式多條染色體並行augustus

mkdir 01-augustsus && cd 01-augustsus
ln ../00-RepeatMask/chi_unmasked.fa.masked genome.fa
seqkit split genome.fa #結果檔案在genome.fa.split
find genome.fa.split/ -type f -name "*.fa" | parallel -j 30 augustus --species=arabidopsis --gff3=on >> temp.gff #並行處理
join_aug_pred.pl < temp.gff  | grep -v '^#' > temp.joined.gff
bedtools sort -i temp.joined.gff > augustsus.gff

AUGUSTUS依賴於已有的模型，而GeneMark-ES/ET則是唯一一款支援無監督訓練模型，之後再識別真核基因組蛋白編碼區的工具。

gmes_petap.pl --ES --sequence genome.fa --cores 50

最後得到的是genemark.gtf，是標準的GTF格式，可以使用Sequence Ontology Project提供的gtf2gff3.pl進行轉換

wget http://genes.mit.edu/burgelab/miso/scripts/gtf2gff3.pl
chmod 755 gtf2gff3.pl
gtf2gff3.pl genemark.gtf | bedtools sort -i - > genemark.gff

不同從頭預測軟體的實際效果可以通過在IGV中載入文章提供的gff檔案和預測後的gff檔案進行比較，一般會存在如下幾個問題：

基因多了，或者少了，也就是假陽性和假陰性現象
UTR區域難以預測，這個比較正常
未正確識別可變剪下位點，導致前後幾個基因識別成一個基因

考慮到轉錄組測序已經非常便宜，可以通過該物種的RNA-seq提供覆蓋度資訊進行預測。

基於轉錄組資料預測

根據已有的模型或者自訓練可以正確預測很大一部分的基因，但如果需要提高預測的正確性，還需要額外的資訊。在過去就需要提供物種本身的cDNA, EST，而現在更多的是基於轉錄組序列進行訓練。儘管RNA-seq資料在基因組上的比對情況能夠推測出內含子位置，根據覆蓋度可以推測出外顯子和非編碼區的邊界，但是僅僅依賴於RNA-seq的覆蓋不能可信地推測出蛋白編碼區(Hoff K.J. Stanke M. 2015).

AUGUSTUS可以利用轉錄組比對資料中的位置資訊來訓練模型，GeneMark-ET可以利用RNA-seq得到的內含子位點資訊自我訓練HMM引數，進行基因預測。BRAKER2將兩者進行整合，使用GeneMark-ET根據RNA-seq無監督訓練模型尋找基因，然後用AUGUSTUS進行模型訓練，最後完成基因預測

BRAKER流程

首先使用hisat2根據遮蔽後的參考序列建立索引，進行比對。

# 專案根目錄
mkdir index
hisat2-build 01-augustus/genome.fa index/chi_masked
hisat2 -p 20 -x index/chi_masked -1 rna-seq/leaf_ox_r1_1.fastq.gz -2 rna-seq/leaf_ox_r1_2.fastq.gz | samtools sort [email protected] 10 > 02-barker/leaf_ox_r1.bam &
isat2 -p 20 -x index/chi_masked -1 rna-seq/ox_flower9_rep1_1.fastq.gz -2 rna-seq/ox_flower9_rep1_2.fastq.gz | samtools sort [email protected] 10 > 02-barker/ox_flower9.bam &
hisat2 -p 20 -x index/chi_masked -1 rna-seq/ox_flower16_rep1_1.fastq.gz -2 rna-seq/ox_flower16_rep1_2.fastq.gz | samtools sort [email protected] 10 > 02-barker/ox_flower16.bam &

然後，以未遮蔽重複序列的參考序列和BAM檔案作為輸入，讓BRAKER2（安裝會稍顯麻煩，因為依賴許多軟體）進行預測。

braker.pl --gff3 --cores 50 --species=carhr --genome=chi_unmasked.fa --bam=02-barker/leaf_ox_r1.bam,02-barker/ox_flower16.bam,02-barker/ox_flower9.bam
# --gff3: 輸出GFF3格式
# --genome: 基因組序列
# --bam: 比對後的BAM檔案，允許多個
# --cores: 處理核心數

最後會得到如下輸出檔案

hintsfile.gff: 從RNA-seq比對結果的BAM檔案中提取，其中內含子用於訓練GeneMark-EX, 使用所有特徵訓練AUGUSTUS
GeneMark-ET/genemark.gtf: GeneMark-EX根據RNA-seq資料訓練後預測的基因
augustus.hints.gff: AUGUSTUS輸出檔案

將augustus.hints.gff3和文章的註釋檔案(carhr38.gtf)比較，見下圖：

comparation

其實不難發現，在不考慮UTR區域情況下，兩者的差別其實更多表現是基因數目上，其實也就是利用轉錄組資料推測結構的問題所在，沒有覆蓋的區域到底是真的沒有基因，還是有基因結構只不過所用組織沒有表達，或者說那個區域其實是假基因？此外，如果基因間隔區域很短，有時候還會錯誤地把兩個不同的基因預測為一個基因。因此，應該注重RNA-seq資料在剪下位點識別和外顯子邊界確定的優勢。

03-同源預測基因結構

同源預測(homology prediction)利用近緣物種已知基因進行序列比對，找到同源序列。然後在同源序列的基礎上，根據基因訊號如剪下訊號、基因起始和終止密碼子對基因結構進行預測，如下示意圖：

同源註釋

相對於從頭預測的“大海撈針”，同源預測相當於先用一塊磁鐵在基因組大海中縮小了可能區域，然後從可能區域中鑑定基因結構。在10年之前，當時RNA-seq還沒有普及, 只有少部分物種才有EST序列和cDNA序列的情況下，這的確是一個比較好的策略，那麼問題來了，現在還需要進行這一步嗎，如果需要是出於那種角度考慮呢?

在同源預測上，目前看到的大部分基因組文章都是基於TBLASTN + GeneWise，這可能是因為大部分基因組文章都是國內做的，這些註釋自然而言用的就是公司的流程，然後目前國內的公司大多數又和某一家公司有一些關係。不過最近的3010水稻泛基因組用的是MAKER, 感謝部分提到這部分工作是由M. Roa(Philippine Genome Center Core Facilities for Bioinformatics, Department of Science)做的，算是一股清流吧。當然我在看Cardamine hirsuta基因組註釋文章，發現它們同源註釋部分用的是GenomeThreader, 該工具在本篇文章成文時的3月之前又更新了。

GeneWise的網站說它目前由Ewan Birney維護，只不過不繼續開發了，因為Guy Slater開發Exonerate解決了GeneWise存在的很多問題，並且速度快了1000倍。考慮到目前只有GeneWise能利用HMM根據蛋白找DNA，而且ENSEMBL的註釋流程也有一些核心模組用到了它，所以作者依舊在緩慢的開發這個工具(自2.4.1已經10多年沒有更新了)，當然這個工具也是非常的慢。儘管這一步不會用到GeneWise作為我們的同源註釋選項，但是我們可以嘗試用GeneWise手工註釋一個基因，主要步驟如下

第一步：使用BLASTX，根據dna序列搜尋到蛋白序列，只需要第一個最佳比對結果
第二步：選擇最佳比對的氨基酸序列
第三步：將dna序列前後延長2kb，與氨基酸序列一併傳入給genewise進行同源預測

提取前5K序列，然後選擇在TAIR上用BLASTX進行比對

seqkit faidx chi_unmasked.fa Chr1:1-5000 > chr1_5k.fa

BLASTX

選擇第一個比對結果中的氨基酸序列，和前5k的DNA序列一併作為GeneWise的輸入

GeneWise2

最後的結果出乎了我的意料

預測結果

讓我們跳過這個尷尬的環節，畢竟很可能是我不太熟練使用工作所致。這裡說點我的看法，除非你真的沒有轉錄組資料，必須要用到同源物種的蛋白進行預測，或者你手動處理幾個基因，否則不建議使用這個工具，因為你可能連安裝都搞不定。

讓我們用GenomeThreader基於上面的DNA序列和氨基酸序列進行同源基因結構預測吧

gth -genomic chr1_5k.fa -protein cer.fa -intermediate -gff3out
# 其中cer.fa就是AT1G02205.2的氨基酸序列

結果一致，並且從RNA-seq的覆蓋情況也符合預期

Chr1    gth exon    1027    1197    Parent=gene1    Chr1    MIPS_CARH_v3.8  exon    1027    1197
Chr1    gth exon    1275    1448    Parent=gene1    Chr1    MIPS_CARH_v3.8  exon    1275    1448
Chr1    gth exon    1541    1662    Parent=gene1    Chr1    MIPS_CARH_v3.8  exon    1555    1662
Chr1    gth exon    1807    2007    Parent=gene1    Chr1    MIPS_CARH_v3.8  exon    1807    2007
Chr1    gth exon    2085    2192    Parent=gene1    Chr1    MIPS_CARH_v3.8  exon    2085    2192
Chr1    gth exon    2294    2669    Parent=gene1    Chr1    MIPS_CARH_v3.8  exon    2294    2669
Chr1    gth exon    3636    3855    Parent=gene1    Chr1    MIPS_CARH_v3.8  exon    3636    3855
Chr1    gth exon    3971    4203    Parent=gene1    Chr1    MIPS_CARH_v3.8  exon    3971    4203
Chr1    gth exon    4325    4548    Parent=gene1    Chr1    MIPS_CARH_v3.8  exon    4325    4548
Chr1    gth exon    4676    4735    Parent=gene1    Chr1    MIPS_CARH_v3.8  exon    4676    4735

全基因組範圍預測流程如下：

準備cDNA和或protein序列：在https://phytozome.jgi.doe.gov/p下載靠譜的物種的蛋白質序列，如 Arabidopsis thaliana, Oryza sativa, Brassica rapa, 查詢文獻尋找目前該物種的已有EST/cDNA序列，或者RNA-seq從頭組裝轉錄組。這裡僅考慮用同源物種的蛋白序列進行比對分析，轉錄組從頭組裝資料用於PASA整體比對到參考基因組和更新已有的基因解僱。

分別測試下不同物種的同源註釋結果

#run seperately
gth -species arabidopsis -translationtable 1 -gff3 -intermediate -protein ~/db/protein_db/Athaliana_167_TAIR10.protein.fa.gz -genomic chi_unmasked.fa -o 03-genomethreader/Athaliana.gff3 &
gth -species arabidopsis -translationtable 1 -gff3 -intermediate -protein ~/db/protein_db/BrapaFPsc_277_v1.3.protein.fa.gz -genomic chi_unmasked.fa -o 03-genomethreader/Brapa.gff3 &
gth -species arabidopsis -translationtable 1 -gff3 -intermediate -protein ~/db/protein_db/Osativa_323_v7.0.protein.fa.gz -genomic chi_unmasked.fa -o 03-genomethreader/Osativa.gff3 &

在定性角度上來看，同源註釋的結果和從頭預測的沒啥差別, 其中B. rapa和A. thaliana和C. hirsuta都屬於十字花科，而O. sativa是禾本科, 所以前兩者預測的效果好。

IGV展示

當然實際的同源註釋流程中不能是單個物種分別預測，應該是將所有的蛋白序列進行合併，然後用BLASTX找到最優的聯配，之後用GenomeThreader進行預測。PASA流程提到的UniRef90作為同源註釋的搜尋資料庫可能是更好的選擇，由於UniRef優先選擇哪些人工審查、註釋質量高、來源於模式動植物的蛋白，所以可靠性相對於直接使用同源物中可能更高。

BLASTX + GenomeThreader的程式碼探索中

04-RNA-seq的兩種使用策略

對於RNA-seq資料，有兩種使用策略，一種是使用HISAT2 + StringTie先比對再組裝, 一種是從頭組裝，然後使用PASA將轉錄本比對到基因組上。

基於HISAT2 + StringTie

首先，使用HISAT2將RNA-seq資料比對到參考基因組, 這一步和之前相似，但是要增加一個引數--dta，使得StingTie能更好的利用雙端資訊

hisat2-build 01-augustus/genome.fa index/chi_masked
hisat2 --dta -p 20 -x index/chi_masked -1 rna-seq/leaf_ox_r1_1.fastq.gz -2 rna-seq/leaf_ox_r1_2.fastq.gz | samtools sort [email protected] 10 > rna-seq/leaf_ox_r1.bam &
hisat2 --dta -p 20 -x index/chi_masked -1 rna-seq/ox_flower9_rep1_1.fastq.gz -2 rna-seq/ox_flower9_rep1_2.fastq.gz | samtools sort [email protected] 10 > rna-seq/ox_flower9.bam &
hisat2 --dta -p 20 -x index/chi_masked -1 rna-seq/ox_flower16_rep1_1.fastq.gz -2 rna-seq/ox_flower16_rep1_2.fastq.gz | samtools sort [email protected] 10 > rna-seq/ox_flower16.bam &
samtools merge [email protected] 10 rna-seq/merged.bam rna-seq/leaf_ox_r1.bam rna-seq/ox_flower9.bam rna-seq/ox_flower16.bam

然後用StringTie進行轉錄本預測

stringtie -p 10 -o rna-seq/merged.gtf rna-seq/merged.bam

對於後續的EvidenceModeler而言，它不需要UTR資訊，只需要編碼區CDS，需要用TransDecoder進行編碼區預測

util/cufflinks_gtf_genome_to_cdna_fasta.pl merged.gtf input/chi_masked.fa > transcripts.fasta
util/cufflinks_gtf_to_alignment_gff3.pl merged.gtf > transcripts.gff3
TransDecoder.LongOrfs -t transcripts.fasta
TransDecoder.Predict -t transcripts.fasta
util/cdna_alignment_orf_to_genome_orf.pl \
     transcripts.fasta.transdecoder.gff3 \
     transcripts.gff3 \
     transcripts.fasta > transcripts.fasta.transdecoder.genome.gff3

最後結果transcripts.fasta.transdecoder.gff3用於提供給EvidenceModeler

基於PASA

在多年以前，那個基因組組裝還沒有白菜價，只有幾個模式物種基因組的時代，對於一個未測序的基因組，研究者如果要研究某一個基因的功能，大多會通過同源物種相似基因設計PCR引物，然後去擴增cDNA. 如果是一個已知基因組的物種，如果要大規模識別基因, 研究者通常會使用EST(expressed sequence tags)序列。

相對於基於演算法的從頭預測，cDNA和EST序列更能夠真實的反應出一個基因的真實結構，如可變剪下、UTR和Poly-A位點。PASA(Progam to Assemble Spliced Alignments)流程最早用於擬南芥基因組註釋，最初的設計是通過將全長(full-length)cDNA和EST比對到參考基因組上，去發現和更新基因組註釋。其中FL-cDNA和EST序列對最後結果的權重不同。

PASA流程示意

這是以前的故事，現在的故事是二代轉錄組以及一些三代轉錄組資料，那麼如何處理這些資料呢？我認為三代轉錄組相對於過去的FL-cDNA，而二代轉錄組資料經過拼接後可以看作是更長的EST序列。由於目前最普及的還是普通的mRNA-seq, 也就只介紹這部分流程。

考慮到我還沒有研究過三代的全長轉錄組，分析過資料，這裡的思考極有可能出錯，後續可能會修改這一部分思考。

轉錄組組裝使用Trinity(conda安裝)

cd rna-seq
Trinity --seqType fq --CPU 50 --max_memory 64G --left leaf_ox_r1_1.fastq.gz,ox_flower16_rep1_1.fastq.gz,ox_flower9_rep1_1.fastq.gz --right leaf_ox_r1_2.fastq.gz,ox_flower16_rep1_2.fastq.gz,ox_flower9_rep1_2.fastq.gz &

PASA是由30多個命令組成的流程，相關命令位於PASApipeline/scripts，為了適應不同的分析，有些引數需要通過修改配置檔案更改,

cp ~/opt/biosoft/PASApipeline/pasa_conf/pasa.alignAssembly.Template.txt alignAssembly.config
# 修改如下內容
DATABASE=database.sqlite
validate_alignments_in_db.dbi:--MIN_PERCENT_ALIGNED=80
validate_alignments_in_db.dbi:--MIN_AVG_PER_ID=80

上述幾行配置檔案表明SQLite3資料庫的名字，設定了scripts/validate_alignments_in_db.dbi的幾個引數, 表示聯配程度和相似程度。後續以Trinity組裝結果和參考基因組作為輸入，執行程式：

~/opt/biosoft/PASApipeline/scripts/Launch_PASA_pipeline.pl -c alignAssembly.config -C -R -g ../chi_unmasked.fa -t ../rna-seq/trinity_out_dir/Trinity.fasta --ALIGNERS blat,gmap

最後結果如下：

database.sqlite.pasa_assemblies_described.txt
database.sqlite.pasa_assemblies.gff3
database.sqlite.pasa_assemblies.gtf
database.sqlite.pasa_assemblies.bed

其中gff3格式用於後續的分析。

目前的一些想法，將從頭組裝的轉錄本比對到參考基因組上很大依賴組裝結果，所以和EST序列和cDNA相比，質量上還有一點差距。

05-整合預測結果

從頭預測，同源註釋和轉錄組整合都會得到一個預測結果，相當於收集了大量證據，下一步就是通過這些證據定義出更加可靠的基因結構，這一步可以通過人工排查，也可以使用EVidenceModeler(EVM). EVM只接受三類輸入檔案：

gene_prediction.gff3: 標準的GFF3格式，必須要有gene, mRNA, exon, CDS這些特徵，用EVidenceModeler-1.1.1/EvmUtils/gff3_gene_prediction_file_validator.pl驗證
protein_alignments.gff3: 標準的GFF3格式，第9列要有ID信和和target資訊, 標明是比對結果
transcript_alignments.gff3:標準的GFF3格式，第9列要有ID信和和target資訊，標明是比對結果

EVM對gene_prediction.gff3有特殊的要求，就是GFF檔案需要反映出一個基因的結構，gene->(mRNA -> (exon->cds(?))(+))(+), 表示一個基因可以有多個mRNA，即基因的可變剪接, 一個mRNA都可以由一個或者多個exon(外顯子), 外顯子可以是非翻譯區(UTR),也可以是編碼區(CDS). 而GlimmerHMM, SNAP等

這三類根據人為經驗來確定其可信度，從直覺上就是用PASA根據mRNA得到的結果高於從頭預測。

第一步：建立權重檔案,第一列是來源型別(ABINITIO_PREDICTION, PROTEIN, TRANSCRIPT), 第二列對應著GFF3檔案的第二列，第三列則是權重.我這裡用了三個來源的資料。

mkdir 05-EVM && cd 05-EVM
#vim weights.txt
ABINITIO_PREDICTION      augustus       4
TRANSCRIPT      assembler-database.sqlite      7
OTHER_PREDICTION  transdecoder  8

我覺得根據基因組引導組裝的ORF的可信度高於組裝後比對，所以得分和PASA差不多一樣高。從頭預測權重一般都是1，但是BRAKER可信度稍微高一點，可以在2~5之間。

第二步：分割原始資料, 用於後續並行. 為了降低記憶體消耗，–segmentsSize設定的大小需要少於1Mb(這裡是100k)， –overlapSize的不能太小，如果數學好，可用設定成基因平均長度加上2個標準差，數學不好，就設定成10K吧

cat transcripts.fasta.transdecoder.genome.gff3 ../braker/carhr/augustus.hints.gff3 > gene_predictions.gff3
ln ../04-align-transcript/database.sqlite.pasa_assemblies.gff3 transcript_alignments.gff3
~/opt/biosoft/EVidenceModeler-1.1.1/EvmUtils/partition_EVM_inputs.pl --genome ../chi_unmasked.fa --gene_predictions gene_predictions.gff3 --transcript_alignments transcript_alignments.gff3 --segmentSize 100000 --overlapSize 10000 --partition_listing partitions_list.out

第三步：建立並行運算命令並且執行

~/opt/biosoft/EVidenceModeler-1.1.1/EvmUtils/write_EVM_commands.pl --genome ../chi_unmasked.fa --weights `pwd`/weights.txt \
      --gene_predictions gene_predictions.gff3 \
      --transcript_alignments transcript_alignments.gff3 \
      --output_file_name evm.out  --partitions partitions_list.out >  commands.list
parallel --jobs 10 < commands.list

第四步：合併並行結果

~/opt/biosoft/EVidenceModeler-1.1.1/EvmUtils/recombine_EVM_partial_outputs.pl --partitions partitions_list.out --output_file_name evm.out

第五步：結果轉換成GFF3

~/opt/biosoft/EVidenceModeler-1.1.1/EvmUtils/convert_EVM_outputs_to_GFF3.pl  --partitions partitions_list.out --output evm.out  --genome ../chi_unmasked.fa
find . -regex ".*evm.out.gff3" -exec cat {} \; | bedtools sort -i - > EVM.all.gff

當前權重設定下，EVM的結果更加嚴格，需要按照實際情況調整，增加其他證據。

06-可選步驟

註釋過濾：對於初步預測得到的基因，還可以稍微優化一下，例如剔除編碼少於50個AA的預測結果，將轉座子單獨放到一個檔案中(軟體有TransposonPSI)。

這裡基於gffread先根據註釋資訊提取所有的CDS序列，過濾出長度不足50AA的序列，基於這些序列過濾原來的的註釋

gffread EVM.all.gff -g input/genome.fa -y tr_cds.fa
bioawk -c fastx '$seq < 50 {print $comment}' tr_cds.fa | cut -d '=' -f 2 > short_aa_gene_list.txt
grep -v -w -f short_aa_gene_list.txt EvM.all.gff > filter.gff

使用PASA更新EVM結果：EVM結果不包括UTR區域和可變剪下的註釋資訊，可以使用PASA進行更新。然而這部分已經無法逃避MySQL, 伺服器上並沒有MySQL的許可權，我需要學習Perl指令碼進行修改。因此基因結構註釋到此先放一放。

07-基因編號

對每個基因實現編號，形如ABCD000010的效果，方便後續分析。如下程式碼是基於EVM.all.gff，使用方法為python gffrename.py EVM_output.gff prefix > renamed.gff.

#!/usr/bin/env python3
import re
import sys

if len(sys.argv) < 3:
    sys.exit()

gff = open(sys.argv[1])
prf = sys.argv[2]

count = 0
mRNA  = 0
cds   = 0
exon  = 0

print("##gff-version 3.2.1")
for line in gff:
    if not line.startswith("\n"):
        records = line.split("\t")
        records[1] = "."
    if re.search(r"\tgene\t", line):
        count = count + 10
        mRNA  = 0
        gene_id = prf + str(count).zfill(6)
        records[8] = "ID={}".format(gene_id)
    elif re.search(r"\tmRNA\t", line):
        cds   = 0
        exon  = 0
        mRNA  = mRNA + 1
        mRNA_id    = gene_id + "." + str(mRNA)
        records[8] = "ID={};Parent={}".format(mRNA_id, gene_id)
    elif re.search(r"\texon\t", line):
        exon     = exon + 1
        exon_id  = mRNA_id + "_exon_" + str(exon)
        records[8] = "ID={};Parent={}".format(exon_id, mRNA_id)
    elif re.search(r"\tCDS\t", line):
        cds     = cds + 1
        cds_id  = mRNA_id + "_cds_" + str(cds)
        records[8] = "ID={};Parent={}".format(cds_id, mRNA_id)
    else:
        continue

    print("\t".join(records))

gff.close()

一些經驗

如果有轉錄組資料，沒必須要使用太多的從頭預測工具，braker2 加 GlimmerHMM可能就夠用了, 更多是使用PASA和StringTie利用好轉錄組資料進行註釋。

基因功能註釋

基因功能的註釋依賴於上一步的基因結構預測，根據預測結果從基因組上提取翻譯後的 蛋白序列 和主流的資料庫進行比對，完成功能註釋。常用資料庫一共有以幾種：

Nr：NCBI官方非冗餘蛋白資料庫，包括PDB, Swiss-Prot, PIR, PRF; 如果要用DNA序列，就是nt庫
Pfam: 蛋白結構域註釋的分類系統
Swiss-Prot: 高質量的蛋白資料庫，蛋白序列得到實驗的驗證
KEGG: 代謝通路註釋資料庫.
GO: 基因本體論註釋資料庫

除了以上幾個比較通用的資料庫外，其實還有很多小眾資料庫，應該根據課題研究和背景進行選擇。注意，資料庫本身並不能進行註釋，你只是通過序列相似性進行搜尋，而返回的結果你稱之為註釋。因此資料庫和搜尋工具要進行區分，所以你需要單獨下載資料庫和搜尋工具，或者是同時下載包含資料庫和搜尋工具的安裝包。

注意，後續分析中一定要保證你的蛋白序列中不能有代表氨基酸字元以外的字元，比如說有些軟體會把最後一個終止密碼子翻譯成”.”或者”*”

BLASTP

這一部分用到的資料庫都是用BLASTP進行檢索，基本都是四步發：下載資料庫，構建BLASTP索引，資料庫檢索，結果整理。其中結果整理需要根據BLASTP的輸出格式調整。

Nr的NCBI收集的最全的蛋白序列資料庫，但是無論是用NCBI的BLAST還是用速度比較快DIAMOND對nr進行搜尋，其實都沒有利用好物種本身的資訊。因此在RefSeq上下載對應物種的蛋白序列, 用BLASTP進行註釋即可。

# download
wget -4  -nd -np -r 1 -A *.faa.gz ftp://ftp.ncbi.nlm.nih.gov/refseq/release/plant/
mkdir -p ~/db/RefSeq
zcat *.gz > ~/db/RefSeq/plant.protein.faa
# build index
~/opt/biosoft/ncbi-blast-2.7.1+/bin/makeblastdb -in plant.protein.faa -dbtype prot -parse_seqids -title RefSeq_plant -out plant
# search
~/opt/biosoft/ncbi-blast-2.7.1+/bin/blastp -query protein.fa -out RefSeq_plant_blastp.xml -db ~/db/RefSeq/uniprot_sprot.fasta -evalue 1e-5 -outfmt 5 -num_threads 50 &

Swiss-Prot裡收集了目前可信度最高的蛋白序列，一共有55w條記錄，資料量比較小，

# download
wget -4 -q ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz
gzip -d uniprot_sprot.fasta.gz
# builid index
~/opt/biosoft/ncbi-blast-2.7.1+/bin/makeblastdb -in uniprot_sprot.fasta -dbtype prot -title swiss_prot -parse_seqids
# search
~/opt/biosoft/ncbi-blast-2.7.1+/bin/blastp -query protein.fa -out swiss_prot.xml -db ~/db/swiss_prot/uniprot_sprot.fasta -evalue 1e-5 -outfmt 5 -num_threads 50 &

InterProScan

下面介紹的工具是InterProScan, 從它的9G的體量就可以感受它的強大之處，一次運行同時實現多個資訊註釋。

InterPro註釋
Pfam資料庫註釋(可以通過hmmscan搜尋pfam資料庫完成)
GO註釋(可以基於NR和Pfam等資料庫，然後BLAST2GO完成,)
Reactome通路註釋，不同於KEGG

命令如下

./interproscan-5.29-68.0/interproscan.sh -appl Pfam  -f TSV -i sample.fa -cpu 50 -b sample -goterms -iprlookup -pa

-appl告訴軟體要執行哪些資料分析，勾選的越多，分析速度越慢，Pfam就行。

KEGG

KEGG資料庫目前本地版收費，線上版收費，所以只能將蛋白序列在KEGG伺服器上執行。因此你需要在http://www.genome.jp/tools/kaas/選擇合適的工具進行後續的分析。我上傳的50M大小蛋白序列，在KEGG伺服器上只需要執行8個小時，也就是晚上提交任務，白天回來幹活。

執行時間

附錄

基因組註釋的常用軟體：

重複區域
- RepeatMasker：識別基因組中的可能重複
- RepeatModeler: 識別新的重複序列
從頭預測
- Augustus
- Fgenesh
同源預測
- GeneWise
- Exonerate
- Trinity
- GenomeThreader
註釋合併
- GLEAN：已經落伍於時代了
- EvidenceModeler：與時俱進
流程
- PASA：真核生物基因的轉錄本可變剪下自動化註釋專案，需要提供物種的EST或RNA-seq資料
- MAKER
- BRAKER1: 使用GeneMark-ET和AUGUSTUS基於RNA-Seq註釋基因結構
- EuGene
視覺化
- IGV
- JBrowse/GBrowse

參考文獻和推薦閱讀：

真核基因組註釋入門: “A beginner’s guide to eukaryotic genome annotation”
二代測序註釋流程:Comparative Gene Finding: “Annotation Pipelines for Next-Generation Sequencing Projects”
基因組轉錄組註釋策略: “Plant genome and transcriptome annotations: from misconceptions to simple solution”
重複序列綜述: “Repetitive DNA and next-generation sequencing: computational challenges and solutions”
《生物資訊學》樊龍江: 第1-5章: 基因預測與功能註釋
《NGS生物資訊分析》陳連福：真核生物基因組基因註釋

環境準備

資料下載

# Cardamine hirsutat基因組資料
mkdir chi_annotation && cd chi_annotation
wget http://chi.mpipz.mpg.de/download/sequences/chi_v1.fa
cat chi_v1.fa | tr 'atcg' 'ATCG' > chi_unmasked.fa
# Cardamine hirsutat轉錄組資料
wget -4 -q -A '*.fastq.gz' -np -nd -r 2 http://chi.mpipz.mpg.de/download/fruit_rnaseq/cardamine_hirsuta/ &
wget -4 -q -A '*.fastq.gz' -np -nd -r 2 http://chi.mpipz.mpg.de/download/leaf_rnaseq/cardamine_hirsuta/ &

軟體安裝

RepeatMasker: 用於註釋基因組的重複區，需要安裝RMBlast, TRF，以及在http://www.girinst.org註冊以下載Repbase

安裝RepeatMasker

cd ~/src
wget http://tandem.bu.edu/trf/downloadstrf409.linux64
mv trf409.linux64 ~/opt/bin/trf
chmod a+x ~/opt/bin/trf
# RMBlast
cd ~/src
wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.6.0/ncbi-blast-2.6.0+-src.tar.gz
wget http://www.repeatmasker.org/isb-2.6.0+-changes-vers2.patch.gz
tar xf ncbi-blast-2.6.0+-src
gunzip isb-2.6.0+-changes-vers2.patch.gz
cd ncbi-blast-2.6.0+-src
patch -p1 < ../isb-2.6.0+-changes-vers2.patch
cd c++
./configure --with-mt --prefix=~/opt/biosoft/rmblast --without-debug && make && make install
# RepeatMasker
cd ~/src
wget http://repeatmasker.org/RepeatMasker-open-4-0-7.tar.gz
tar xf RepeatMasker-open-4-0-7.tar.gz
mv RepeatMasker ~/opt/biosoft/
cd ~/opt/biosoft/RepeatMasker
## 解壓repbase資料到Libraries下
## 配置RepatMasker
perl ./configure

在上面的基礎上安裝RepeatModel

# RECON
cd ~/src
wget -4 http://repeatmasker.org/RepeatModeler/RECON-1.08.tar.gz
tar xf RECON-1.08.tar.gz
cd RECON-1.08/src
make && make install
cd ~/src
mv RECON-1.08 ~/opt/biosoft
# nesg
cd ~/src
mkdir nesg && cd nesg
wget -4 ftp://ftp.ncbi.nih.gov/pub/seg/nseg/*
make
mv nmerge nseg ~/opt/bin/
# RepeatScout
http://www.repeatmasker.org/RepeatScout-1.0.5.tar.gz
# RepeatModel
wget -4 http://repeatmasker.org/RepeatModeler/RepeatModeler-open-1.0.11.tar.gz
tar xf RepeatModeler-open-1.0.11.tar.gz
mv RepeatModeler-open-1.0.11 ~/opt/biosoft/
cd ~/opt/biosoft/RepeatModeler-open-1.0.11
# 配置
perl ./configure
export PATH=~/opt/biosoft/maker:$PATH

BLAST，BLAST有兩個版本可供選擇, WuBLAST或者NCBI-BLAST，我個人傾向於NCBI-BLAST，並且推薦使用編譯後二進位制版本，因為編譯實在是太花時間了

cd ~/src
wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.7.1+-x64-linux.tar.gz
tar xf ncbi-blast-2.7.1+-x64-linux.tar.gz -C ~/opt/biosoft
# 環境變數
export PATH=~/opt/biosoft/ncbi-blast-2.7.1+/bin:$PATH
# 用於後續的BRAKER2
conda create -n annotation blast=2.2.31

AUGUSTUS: 可以說是最好的預測軟體，使用conda安裝

source activate annotation
conda install augustus=3.3

GeneMark-ES/ET則是唯一一款支援無監督訓練模型, 軟體下載需要登記

cd ~/src
wget http://topaz.gatech.edu/GeneMark/tmp/GMtool_Qg87n/gm_et_linux_64.tar.gz
tar xf gm_et_linux_64.tar.gz
mv gm_et_linux_64/gmes_petap/ ~/opt/biosoft
wget http://topaz.gatech.edu/GeneMark/tmp/GMtool_Qg87n/gm_key_64.gz
gzip -dc gm_key_64.gz > ~/.gm_key
cpan YAML Hash::Merge Logger::Simple Parallel::ForkManager
echo "export PATH=$PATH:~/opt/biosoft/gmes_petap/" >> ~/.bashrc

GlimmerHMM:

cd ~/src
wget -4 ftp://ccb.jhu.edu/pub/software/glimmerhmm/GlimmerHMM-3.0.4.tar.gz
tar xf GlimmerHMM-3.0.4.tar.gz -C ~/opt/biosoft

SNAP: 基因從頭預測工具，在處理含有長內含子上的基因組上表現欠佳

# 安裝
cd ~/src
git clone https://github.com/KorfLab/SNAP.git
cd SNP
make
cd ..
mv SNAP ~/opt/biosoft
# 環境變數
export Zoe=~/opt/biosoft/SNAP/Zoe
export PATH=~/opt/biosoft/SNAP:$PATH

Exnerate 2.2: 配對序列比對工具，提供二進位制版本, 功能類似於GeneWise，能夠將cDNA或蛋白以gao align的方式和基因組序列聯配。

cd ~/src
wget http://ftp.ebi.ac.uk/pub/software/vertebrategenomics/exonerate/exonerate-2.2.0-x86_64.tar.gz
tar xf exonerate-2.2.0-x86_64.tar.gz
mv exonerate-2.2.0-x86_64 ~/opt/biosoft/exonerate-2.2.0
# .bashrc新增環境變數
export PATH=~/opt/biosoft/exonerate-2.2.0:$PATH
# 或
conda install -c bioconda exonerate

GenomeThreader 1.70: 同源預測軟體，1.7.0版本更新於2018年2月

wget -4 http://genomethreader.org/distributions/gth-1.7.0-Linux_x86_64-64bit.tar.gz
tar xf gth-1.7.0-Linux_x86_64-64bit.tar.gz -C ~/opt/biosoft
# 修改.bashrc增加如下行

              
           
              
              
            
            相關推薦
			   
            
            
            
 

    

    
    如何對基因組序列進行註釋
      
							
							
							基因組組裝完成後，或者是完成了草圖，就不可避免遇到一個問題，需要對基因組序列進行註釋。註釋之前首先得構建基因模型，有三種策略：


從頭註釋(de novo prediction)：通過已有的概率模型來預測基因結構，在預測剪下位點和UTR區準確性較低
同源預測( 

  
 

    

    
    @Autowired的使用：推薦對建構函式進行註釋
      
                轉自：http://www.cnblogs.com/acm-bingzi/p/springAutowired.html



　　在編寫程式碼的時候，使用@Autowired註解是，發現IDE報的一個警告，如下：



　　Spring Team recommends "Al 

  
 

    

    
    @Autowired的使用--Spring規範解釋，推薦對建構函式進行註釋
      
							
							
							一
在編寫程式碼的時候，使用@Autowired註解時，發現IDE報的一個警告，如下：
Spring Team recommends "Always use constructor based dependency injection in your beans 

  
 

    

    
    如何對加密pdf進行註釋
      
                非常簡單又有效的註釋加密pdf的方法，只需兩個步驟，是全網目前最簡單的方法，不需要下載任何工具。

如果需要教程請前往：

這裡這裡: https://pan.baidu.com/s/1iVWXkorrMU0tWuaslXqugQ 

pwd是: yqhk

如下載有任何問題 

  
 

    

    
    鏈式基數排序_對整數序列進行排序,低關鍵字優先的基數排序演算法很有創新。
      package cn.itcast.sort;
import java.util.ArrayList;
import java.util.List;
/**
 * 
 * 請對整數序列進行排序。
         隨機產生1000個整數，其中整數的範圍0~9999
         可以用十進位制的每個位為關 

  
 

    

    
    Autowired的使用：推薦對建構函式進行註釋及其他兩種方式
      
                

在編寫程式碼的時候，使用@Autowired註解是，發現IDE報的一個警告，如下：



　　Spring Team recommends "Always use constructor based dependency injection in your beans.  

  
 

    

    
    python學習——讀取染色體長度（七：for循環對染色體序列進行反向互補）
      導入模塊   int   終端   染色體   文件名   循環   open   sys.argv   pan   對fasta文件genome_test.fa中的染色體序列進行反向互補，並輸出到文件genome_test_RC.fa
genome_test.fa
>chr1ATATATATAT& 

  
 

    

    
    @Autowired的使用：推薦對構造函數進行註釋
      you   編寫   pro   user   意思   point   動態   RR   src   


　　在編寫代碼的時候，使用@Autowired註解是，發現IDE報的一個警告，如下：

　　Spring Team recommends "Always use constructor based  

  
 

    

    
    C# 使用 protobuf 進行對象序列化與反序列化
      member   開源項目   serial   all   序列化與反序列化   ace   ogl   serialize   dll   protobuf 是 google的一個開源項目，可用於以下兩種用途：
（1）數據的存儲（序列化和反序列化），類似於xml、json等；
（2）制作網絡通信協議。 

  
 

    

    
    C#對Json資料進行序列化
       
 
 json格式：我們常見的josn格式資料字串有一般都是一對大括號({}),或者兩對大括號的。下面就這2種常見的json格式的資料介紹一些對json格式的操作。 
 先了解下，我說的2種json格式： 
 第一種格式有一對大括號的： 
  
 第二種有兩對大括號的： 
 {

    \"code\" 

  
 

    

    
    Python中使用numpy對序列進行離散傅立葉變換DFT
       
 
 看了大佬對DFT的介紹後感覺離散傅立葉變換對序列訊號的處理還是很有用的， 
 總結下來就是DFT可以增加有限長序列的長度來提高物理解析度。 
 自己用python中的numpy庫實現了一下： 
 其中繪相簿的使用請參考：Python繪圖 
 將有效長度為4的單位序列，變換為長度16的DFT譜線。 
 

  
 

    

    
    調用FFmpeg SDK對YUV視頻序列進行編碼
      blog   -a   51cto   fcc   ext   pro   分享   fab   img   
由於作者不習慣該編輯器，只是將本文的截圖貼了出來，詳文見：https://www.yuque.com/docs/share/e2ff9da8-678d-49c6-ac9a-ca0f654d3f73調 

  
 

    

    
    給出一個含有n個數字的序列a1,a2,a3,...an，可以進行以下操作： 一次操作定義為對這個序列的每個數字進行以下兩種改變之一： 1.ai ÷ 2 2.ai × 3 每一次的操作中，必須保證至少有
      
							
							
							JAVA 
給出一個含有n個數字的序列a1,a2,a3,…an，可以進行以下操作： 
一次操作定義為對這個序列的每個數字進行以下兩種改變之一： 
1.ai ÷ 2 
2.ai × 3 
每一次的操作中，必須保證至少有一個數字是第1種改變；並且經過每次操作後，每一 

  
 

    

    
    使用Linq中的Distinct方法對序列進行去重操作
      
                
使用Linq提供的擴充套件方法Distinct可以去除序列中的重複元素。
該方法具有以下兩種過載形式：
（1）public static IEnumerable<TSource> Distinct<TSource>(this IEnumerable& 

  
 

    

    
    如何：對 JSON 資料進行序列化和反序列化
      
                JSON（JavaScript 物件符號）是一種高效的資料編碼格式，可用於在客戶端瀏覽器和支援 AJAX 的 Web 服務之間快速交換少量資料。 
本主題演示如何使用 DataContractJsonSerializer 將 .NET 型別物件序列化為 JSON 編碼資料，然 

  
 

    

    
    序列化工具類（{對實體Bean進行序列化操作.}，{將字節數組反序列化為實體Bean.}）
      fin   pub   字節數   字節   工具類   ktr   null   pan   port   
package com.dsj.gdbd.utils.serialize;

import java.io.ByteArrayInputStream;
import java.io.Byte 

  
 

    

    
    序列化工具類（{對實體Bean進行序列化操作.}，{將位元組陣列反序列化為實體Bean.}）
       
 package com.dsj.gdbd.utils.serialize;

import java.io.ByteArrayInputStream;
import java.io.ByteArrayOutputStream;
import java.io.IOException;
impor 

  
 

    

    
    SpringMVC註解@Autowired和@Qualifier  自動注入[根據型別注入]     @Autowired 可以對成員變數、方法以及建構函式進行註釋，     @Qualifier 的
      
                
@Autowired和@Qualifier  自動注入[根據型別注入]  
  @Autowired 可以對成員變數、方法以及建構函式進行註釋，  
  @Qualifier 的標註物件是成員變數、方法入參、建構函式入參。  
  ps：兩者結合使用相當於@Resource 

  
 

    

    
    下載基因組註釋gtf檔案和下載參考基因組序列
      
                

下載GTF註釋檔案



對NCBI：

wget ftp://ftp.ncbi.nlm.nih.gov/genomes/Homo_sapiens/GFF/ref_GRCh38.p7_top_level.gff3.gz          ## hg38
wget ftp: 

  
 

    

    
    【FFMpeg視訊開發與應用基礎】二、呼叫FFmpeg SDK對YUV視訊序列進行編碼
      
							
							
							《FFMpeg視訊開發與應用基礎——使用FFMpeg工具與SDK》視訊教程已經在“CSDN學院”上線，視訊中包含了從0開始逐行程式碼實現FFMpeg視訊開發的過程，歡迎觀看！連結地址：FFMpeg視訊開發與應用基礎——使用FFMpeg工具與SDK