人類參考基因組

人類參考基因組

一、人類參考基因組的來源

1、人類基因組計劃

1）2001年草圖，繪制人類基因組圖譜

2、資料庫的名稱

1）UCSC：hg19，hg38

2）NCBI：GRCH19，GRCH38

二、如何下載下傳參考基因組

在 linux 中下載下傳參考序列資料庫：

1. hg38：wget http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz
2. hg19：wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz 
# 下載下傳時間會比較久，建議網速不好時候，用其他方法。
           

三、參考基因組的後續處理

1、對下載下傳的基因組進行初步整理：

tar zvfx chromFa.tar.gz
# 與 chrNo.fa 無關的檔案删除掉。（chrNo 指 chr1、chr2......chrX、chrY、chrM）

cat  *.fa  >>  genome.fa
rm  chr*.fa    
# 即可獲得人類參考基因組序列的所有染色體的彙總檔案：genome.fa
           

2、用 bwa 軟體，對 genome.fa 建立索引檔案：

[bwa path]  index  genome.fa
# 建構索引後，會生成檔案：hg19.fa.amb、hg19.fa.ann、hg19.fa.bwt、hg19.fa.pac 和 hg19.fa.sa
           

四、參考基因組的資訊統計

1、染色體的長度統計，用 python 畫個條形圖。

答：請參見以下程式1。

2、染色體的 GC 含量，用 python 畫個條形圖。

答：請參見以下程式2。

五、.fai 的深入研究

1、如何得到 .fai 檔案 ?

答：用 samtools 軟體的 faidx 參數，可對參考基因組 .fa 檔案進行索引，生成 .fai 檔案。

例如：[samtools path] faidx input.fa

2、.fai 檔案總共有 5 列資訊：

1）序列所在染色體号

2）序列的長度

3）該序列的第一個堿基，在檔案中的偏移位置，從 0 開始計數，first.base_offset

4）每行有多少個堿基， base_number_each_line

5）每行的位元組長度， byte_each_line

例如：chr1 12345678 68 70 71

示例：

比如，2号染色體的第一個堿基的偏移量 = 1号染色體的長度/每行的堿基數 X 每行的位元組數 + （>chr1\n 和 >chr2\n）的位元組數。

3、如何利用 fai 檔案快速得到某個 region 的參考序列 ?

答：用 samtools 軟體的 faidx 參數，可快得到某個 region 的參考序列。

例如，如果要獲得參考序列 hg19 的1号染色體上的 2377-2399 位置的序列：samtools faidx hg19.fa ’ chr1: 2377-2399 ’

4、給定一個 region，自己用 python 實作參考序列的快速查找 ？

答：利用 seek 函數，用 python 實作參考序列的快速查找。

請參見以下程式3。

程式1：

#!python3

# 讀取genome.fa，進行染色體的長度統計，并用條形圖呈現出來。

import sys
import matplotlib.pyplot as plt
import numpy as np

if len(sys.argv) != 2:
        print('Usage: python3 %s <genome.fa> ' % sys.argv[0])
        sys.exit()

genome_fa = sys.argv[1]
fasta_open = open(genome_fa, "r")
l_hash = {}

# 染色體的長度統計，将染色體号和對應的長度，儲存在字典中。
for line in fasta_open:
    line = line.strip()
    if line.startswith('>'):
        chr_key = line.lstrip('>')      
        l_hash[chr_key] = 0
    else:
        l_hash[chr_key] += len(line)        
fasta_open.close()

chrnum = list(l_hash.keys())
length = list(l_hash.values())

fig1 = plt.figure(num=1, figsize=(8,6))
plt.title('The length of each chromosome', fontweight='bold', fontsize=12)
plt.xlabel('Chromosome', fontsize=9)
plt.ylabel('Number of Base', fontsize=9)
plt.yscale('linear')
plt.xticks(rotation=45, fontsize=7)
plt.bar(chrnum, length)
plt.show()
fig1.savefig('Length_of_Chromosome.png')
           

程式2：

#!python3

# 讀取 genome.fa，進行染色體的GC總含量的統計，并用條形圖呈現出來。

import sys
import matplotlib.pyplot as plt

if len(sys.argv) != 2:
        print('Usage: python3 %s <genome.fa> ' % sys.argv[0])
        sys.exit()

genome_fa = sys.argv[1]
fasta_open = open(genome_fa, "r")

b_hash = {}

# 染色體号為 key，GC 含量為 value，将統計數值儲存在字典中。
for line in fasta_open:
    line = line.strip()
    if line.startswith('>'):
        chr_key = line.lstrip('>')
        b_hash[chr_key] = 0
    else:
        for base in line:
            if base == 'C' or base == 'G':
                b_hash[chr_key] += 1
fasta_open.close()

chrnum = list(b_hash.keys())
gc = list(b_hash.values())

# 繪制條形圖。
fig1 = plt.figure(num=1, figsize=(8,6))
plt.title('The GC contents of each chromosome', fontweight='bold', fontsize=12)
plt.xlabel('Chromosome', fontsize=9)
plt.ylabel('GC Contents', fontsize=9)
plt.yscale('linear')
plt.xticks(rotation=45, fontsize=7)
plt.bar(chrnum, gc)
plt.show()
fig1.savefig('GC_Content_of_Chromosome.png')
           

程式3：

# 實作同 samtools faidx 一樣功能的程式。 

import sys

def Read_Base(filefa, offset, pos2, pos1):
    seek = offset + pos1 - 1    # 該染色體第一個堿基偏移位置 + 序列的第一個堿基在染色體上的位置（-1是為了包括序列的第一個堿基，否則，會從序列的第二個堿基開始）
    pos = pos2 - pos1 + 1   # 序列的長度。
    with open(filefa, 'r') as fin:  # 讀 genome.fa，找序列的堿基段。 
        frontline = pos1 // 50  # 因為 genome.fa 是 50 個堿基為一行，看這條序列前面的堿基跨越了多少行。
        remainline = pos1 % 50  # 看序列的第一個堿基在genome.fa中，是不是剛好在 \n 上（\n算是一個位元組）
        if remainline == 0: seek = seek + frontline - 2  # 如果在 \n 上，則指針要後移2個位置。
        else: seek = seek - 1   # 如果不在 \n 上，指針隻需要後移1個位置。

        seekinitial = seek  # 儲存此時指針的位置，以防下面不滿足規定時，需要回到這個位置上。
        fin.seek(seek, 0)   # 從頭開始，看 seek 位置的堿基。
        text = fin.read(pos + 1)    
        flag = True
        while flag: # 循環程式目的：加上這條序列橫跨的行，即\n的位元組數。
            if '\n' in text:    # 如果該序列包括\n的話。
                seek = seekinitial  # 從原來的指針位置開始讀。              
                count_n = text.count('\n')  # 計算 \n 的數量。
                seek = seek + count_n   # 原來seek需要加上 \n 的位元組數。
                fin.seek(seek, 0)
                text = fin.read(pos + 1)
                if text.count('\n') <= count_n: # 如果seek加上原來\n的位元組數後，又遇到\n，則再循環確定新的\n囊括進來。
                    flag = False
            else:   # 如果該序列不包括\n，則不用加\n的位元組數。
                flag = False
        
        text = text.replace('\n','')  # 去掉\n，隻留下堿基。      
        lenbase = len(text)
        baseline = lenbase // 60   # 模仿 samtools faidx 輸出，每行為60個堿基。便于核對。
        for l in range(baseline+1):
            for t in text[0+60*l:60+60*l]:  # 讀取60個堿基，輸出為同一行。
                print(t, end='')
            print('\n')        
    return 

if __name__ == '__main__':
    if len(sys.argv) != 6:
        print('Usage: python3 run.py <genome.fa> <genome.fa.fai> <which chr> <start pos> <end pos>')
        sys.exit()

    chrnum = sys.argv[3]   # 讀取哪一條染色體的序列。
    num = chrnum.lstrip('chr')  

    # 如果是 chrX chrY chrM，轉化為對應的數字。
    if num == 'X':
        num = 23
    elif num == 'Y':
        num = 24
    elif num == 'M':
        num = 25
    num = int(num)
    pos1 = int(sys.argv[4])
    pos2 = int(sys.argv[5])
    print('>{0}:{1}-{2}'.format(chrnum, pos1, pos2))    # 輸出類似 >chrY:1000-1050 的資訊。

    with open(sys.argv[2], 'r') as f:  # 讀取 genome.fa.fai 檔案。
        for i in range(num):   # 看看是哪一條染色體，則讀取對應行的資訊。
            line = f.readline() # 擷取想要讀取的那一條染色體的資訊。
        linelst = line.strip().split('\t')
        offset = int(linelst[2])
        Read_Base(sys.argv[1], offset, pos2, pos1)