簡介

通常使用bwa做mapping後會獲得sam檔案，而sam檔案包含2個重要的欄位：該序列mapping上的染色體id和位置(比如第2列(chr5)和第3列(36345037))

KMER_44   0     chr5    36345037        37      7M1D24M *       0       0       CTGATGCAAAAAAAAAAAAGCTTTTTTGAAG II*II?IIIIIIIIII8;[email protected];I4E+2 XT:A:U  NM:i:1  X0:i:1  X1:i:0  XM:i:0  XO:i:1  XG:i:1  MD:Z:7^A24
KMER_45   0     chr5    142493459       37      13M2D18M        *       0    
 

若想通過這2個欄位去檢視該位置上包含哪些已經存在的基因，怎麼做呢？
python提供了一個神奇的包：pyensembl
官網demo:

from pyensembl import EnsemblRelease

# release 77 uses human reference genome GRCh38
data = EnsemblRelease(77)

# will return ['HLA-A'] contig就是chr6的6(染色體id)
gene_names = data.gene_names_at_locus(contig=6, position=29945884)

# get all exons associated with HLA-A
 

exon_ids = data.exon_ids_of_gene_name('HLA-A')

安裝步驟

1.先安裝pyensembl

pip install pyensembl

2.執行的時候，提示你要下載一下基因資料庫，執行下面的命令即可

pyensembl install --release 77 --species homo_sapiens

但是不知為何，windowns總是出現如下錯誤:

urlopen error ftp error: error_perm('550 Failed to change directory.',)>

總共需要下載4個壓縮包

Homo_sapiens.GRCh38.77.gtf.gz
Homo_sapiens.GRCh38.cdna.all.fa.gz
Homo_sapiens.GRCh38.ncrna.fa.gz
Homo_sapiens.GRCh38.pep.all.fa.gz

python-ftp-bug無解，只能手動一個一個下載，如下提供一個csdn下載地址:
https://download.csdn.net/download/jiangpeng59/10729670

若是windowns需解壓到如下目錄，注意路徑和自己的使用者名稱對應即可(PJ.Javis)

C:\Users\PJ.Javis\AppData\Local\pyensembl\GRCh38\ensembl77\pyensembl\GRCh38\ensembl77\Cache

然後在此執行上面的命令，但又會提示

ModuleNotFoundError: No module named 'resource'

明明安裝了resource但就是找不到，貌似是win的BUG，如下提供一個暴力的方法，修改檔案"anaconda3\lib\site-packages\gtfparse\util.py"，註釋掉相應的內容- -，最後可正常執行命令。

from __future__ import print_function, division, absolute_import
import sys
#import resource

def memory_usage():
    """
    Returns number of megabytes of memory currently being used by this Python
    process
    """
    # resources = resource.getrusage(resource.RUSAGE_SELF)
    # if sys.platform == 'darwin':
    #     resident_bytes = resources.ru_maxrss
    #     resident_kilobytes = resident_bytes / 1024
    # else:
    #     resident_kilobytes = resources.ru_maxrss
    # return resident_kilobytes / 1024
    return 1024