標籤：

【翻譯】管理高質量的參考基因組

04-18

本秘笈闡述一些管理參考基因組的一般技巧。作為一個直觀的例子，這裡將研究GC含量（基因組中含有鳥嘌呤Guanine-胞嘧啶Cytosine的比例）。參考基因組一般採用FASTA的文件格式。

準備好

不同基因組規模的大小有很大不同，從病毒如艾滋病毒HIV（它有9.7kbp），到細菌如大腸桿菌E.coli，到原生動物如惡性瘧原蟲Plasmodium falciparum（導致瘧疾的最重要的寄生物種）有14條染色體、線粒體和/*譯者註：一個高度退化的質體樣細胞器 */ apicoplast，到果蠅有三個常染色體、線粒體和X/Y性染色體，到人類有3Gbps分布於22個常染色體，X/Y染色體和線粒體中，最後直到巴黎粳稻Paris japonica，一種具有150 Gbp基因組的植物。 一路走來，讀者可以看到有不同的染色體倍性和不同的性染色體的組織機構。

小技巧

正如你看到的，不同物種有不同大小的基因組。這種不同可以有好幾個數量級。這對用戶的編程風格有顯著的影響。用大的基因組工作需要用戶更小心地處理內存的使用，不幸的是，有時間效率的編程技術對更大的基因組也會有更大的收益（因為你有更多的數據需要分析），這些是相互矛盾的需求。一般的規則是用戶必須在面對大的基因組時做得效率更高（包括速度和內存）。

為了使這個秘笈減輕負擔，我們將使用來自惡性瘧原蟲的小型真核基因組。該基因組仍然具有較大基因組的許多典型特徵（例如，多個染色體），所以這是複雜性和大小之間的一個很好的折衷。請注意，對於惡性瘧原蟲大小的基因組，通過將全基因組載入到內存中將可能執行許多操作。然而，我們選擇了一種可以與更大的基因組（例如哺乳動物）一起使用的編程風格，以便讀者能以更一般的方式使用此秘笈，但讀者可以隨意使用更多的內存密集型方法來處理像這樣的小型基因組。

我們將使用第1章中安裝的Biopython，第一章，Python和周邊的軟體生態，像往常一樣，本秘笈可以在IPython Notebook的02_Genomes/Reference_Genome.ipynb 找到。

如果讀者沒有使用notebook可以從我們的數據集中下載P. falciparum的基因組，網頁在github.com/tiagoantao/b (文件 pfalciparum.fasta)

如何做

參看以下步驟：

1.這裡先檢查參考基因組FASTA文件中所有序列的描述：

from Bio import SeqIOgenome_name = PlasmoDB-9.3_Pfalciparum3D7_Genome.fastarecs = SeqIO.parse(genome_name, fasta)chroms = {}for rec in recs: print(rec.description)

  • 這個代碼看起來和前一章類似；看一下部分的輸出：

  • 不同基因組的參考序列將具有不同的描述行，但是它們一般會在那裡包括重要的信息。在這個例子中，讀者可以看到有chromosomes, mitochondria和apicoplas。這裡也可以查看染色體大小，但是我們將從序列長度中取值。

2. 解析這個描述行來提取染色體的序號。我們將從序列中取得染色體的大小，然後按窗口計算染色體上的GC含量：

from Bio import SeqUtilschrom_sizes = {}chrom_GC = {}recs = SeqIO.parse(genome_name, fasta)block_size = 50000min_GC = 100.0max_GC = 0.0for rec in recs: if rec.description.find(SO=chromosome) == -1: continue chrom = int(rec.description.split(_)[1]) chrom_GC[chrom] = [] size = len(rec.seq) chrom_sizes[chrom] = size num_blocks = size // block_size + 1 for block in range(num_blocks): start = block_size * block if block == num_blocks - 1: end = size else: end = block_size + start + 1 block_seq = rec.seq[start:end] block_GC = SeqUtils.GC(block_seq) if block_GC < min_GC: min_GC = block_GC if block_GC > max_GC: max_GC = block_GC chrom_GC[chrom].append(block_GC)print(min_GC, max_GC)

  • 這裡對所有染色體進行窗口分析，與讀者在前一章中已經看到的類似。從定義一個50kbp的大小的窗口開始，這適用於惡性瘧原蟲（這是可以隨意改變其大小的），但是讀者可以考慮對數量級有所不同的染色體的基因組採用其他值。
  • 請注意，這裡正在重新讀取文件。對於這麼小的基因組，用戶可以（在第一步）靈活地完成在內存中載入整個基因組。不管怎樣，讀者可以隨意對小基因組嘗試這種編程風格- 它更快！但是，對於較大的基因組，這個代碼更為廣泛適用。
  • 請注意，在for循環中，通過解析描述中的SO條目來忽略線粒體和apicoplast。 chrom_sizes字典將存儲染色體的大小。
  • chrom_GC字典是我們最感興趣的數據結構，它是一個列表，每項對應一個50 kbp窗口的GC含量比例。因此，對於大小為640,851bp的染色體1，將有14個數據項，因為該染色體大小具有14個50kbp的區塊。

小技巧

請注意惡性瘧原蟲基因組有兩個不同尋常的特徵：基因組非常富含AT，即缺乏GC。所以，讀者會得到非常低的數字。此外，染色體是根據大小排列的（因為它很常見），但是從最小的開始。通常的慣例是以最大值開始（例如，像人類中的基因組）。

3. 現在實現繪製基因組GC分布的圖，這裡將使用GC含量的藍色陰影。但是，對於高異常值，將使用紅色陰影。 對於低異常值，將使用黃色陰影：

from __future__ import divisionfrom reportlab.lib import colorsfrom reportlab.lib.units import cmfrom Bio.Graphics import BasicChromosome

  • 這裡將用到浮點除法，這需要導入Biopython和reportlab庫

小技巧

在Python成為如此時髦的語言之前，Biopython代碼已經發展了很多年。過去，可用的庫非常有限。 reportlab的用法大多被視為遺留問題。 作者建議讀者從中學到可以與Biopython一起使用的東西即可。 如果讀者打算學習Python的現代繪圖庫，可能需要考慮matplotlib，Bokeh或Python版本的ggplot（或其他可視化替代方法，如Mayavi，VTK或甚至Blender的API）。

chroms = list(chrom_sizes.keys())chroms.sort()biggest_chrom = max(chrom_sizes.values())my_genome = BasicChromosome.Organism(output_format="png")my_genome.page_size = (29.7*cm, 21*cm) # checktelomere_length = 10bottom_GC = 17.5top_GC = 22.0for chrom in chroms: chrom_size = chrom_sizes[chrom] chrom_representation = BasicChromosome.Chromosome(Cr %d % chrom) chrom_representation.scale_num = biggest_chrom tel = BasicChromosome.TelomereSegment() tel.scale = telomere_length chrom_representation.add(tel) num_blocks = len(chrom_GC[chrom]) for block, gc in enumerate(chrom_GC[chrom]): my_GC = chrom_GC[chrom][block] body = BasicChromosome.ChromosomeSegment() if my_GC > top_GC: body.fill_color = colors.Color(1, 0, 0) elif my_GC < bottom_GC: body.fill_color = colors.Color(1, 1, 0) else: my_color = (my_GC - bottom_GC) / (top_GC - bottom_GC) body.fill_color = colors.Color(my_color, my_color, 1) if block < num_blocks - 1: body.scale = block_size else: body.scale = chrom_size % block_size chrom_representation.add(body) tel = BasicChromosome.TelomereSegment(inverted=True) tel.scale = telomere_length chrom_representation.add(tel) my_genome.add(chrom_representation)my_genome.draw("falciparum.png", "Plasmodium falciparum")

  • 第一行用keys方法的返回值轉換為列表。這在Python 2中是多餘的，但在Python 3中不是這樣，其中keys方法具有特定的返回類型：dict_keys。
  • 這裡將按順序繪製染色體（因此排序）。所以需要最大染色體（瘧原蟲中14個）的大小，以確保染色體的大小以正確的比例（最大染色體變數）輸出。
  • 然後創建一個A4大小的PNG表示輸出生物體。請注意，這裡將繪製10 bp的非常小的端粒。這會產生一個矩形的染色體。讀者可以設置更大的端粒，用一個圓角表示（或者你可以對了解的物種的端粒的正確大小有更好的主意）。
  • 這裡聲明GC含量低於17.5％或高於22.0％的任何含量都將被視為異常值。請記住，對於大多數其他物種，這將會更高。
  • 然後，輸出這些染色體的恰當值。它們以端粒為邊界，由50 kbp染色體片段組成（最後一個片段與剩餘片段一起，大小一樣）。基於兩個異常值之間的線性歸一化值，每個區段著相應的藍色，並包含紅色或綠色分量。如果最後一條染色體存在，每條染色體片段可以是50kbp或可能更小。輸出如下圖所示

圖1：惡性瘧原蟲的14條染色體以GC含量進行顏色編碼（紅色大於22％，黃色小於17％，藍色表示兩個數字之間的線性梯度）

4. 最後，如果讀者在IPython NoteBook中（即不是在IPython之外）。你可以內嵌輸出這個圖：

from IPython.core.display import ImageImage("falciparum.png")

附加信息

  • 惡性瘧原蟲是真核生物的一個合理的例子，它具有小的基因組，使讀者可以進行小數據練習，並具有足夠的特徵，使其對大多數真核生物仍然有用。 當然，沒有性染色體（如人類中的X / Y），但這些應該很容易處理，因為參考基因組不會涉及染色體倍性問題。
  • 惡性瘧原蟲確實有線粒體，但由於空間問題，我們不會在這裡處理。 Biopython具有列印輸出環狀基因組的功能，讀者也可以對細菌使用它們。 關於細菌和病毒，這些基因組更容易處理，因為它們的大小非常小。

參考資料

  • 讀者可以在Ensembl網站ensembl.org/info/data/f找到許多模式生物的參考基因組。
  • 通常，NCBI還提供了大量的基因組列表，ncbi.nlm.nih.gov/genome。
  • 有很多網站專門針對單個生物體（或一組相關的生物體）：除了在PlasmoDB（plasmodb.org/plasmo/），這個可以下載了惡性瘧原蟲基因組的資料庫，還可以發現VectorBase（vectorbase.org/）疾病媒介資料庫應用於下一個秘笈。 Flybase（flybase.org/）資料庫對果蠅也值得一提，但讀者不要忘記搜索自己感興趣的有機體。

推薦閱讀：

TAG:Python入門 | 生物信息學 | 繪圖 |