Pengantar Komprehensif Untuk Genom Anda Dengan SciPy Stack

Bioinformatika adalah bidang interdisipliner yang mengembangkan metode dan perangkat lunak untuk menganalisis dan memahami data biologis.

Lebih sederhananya, Anda bisa menganggapnya sebagai ilmu data untuk biologi.

Di antara banyak jenis data biologis, data genomik adalah salah satu yang paling banyak dianalisis. Terutama dengan kemajuan pesat teknologi sekuensing DNA (NGS) generasi berikutnya, volume data genomik telah tumbuh secara eksponensial. Menurut Stephens, Zachary D et al., akuisisi data genomik berada pada skala exabyte per tahun.

Dalam posting ini, saya mendemonstrasikan contoh menganalisis file GFF3 untuk genom manusia dengan SciPy Stack. Format Fitur Generik Versi 3 (GFF3) adalah format file teks standar saat ini untuk menyimpan fitur genomik. Secara khusus, dalam posting ini Anda akan mempelajari cara menggunakan tumpukan SciPy untuk menjawab pertanyaan berikut tentang genom manusia:

1. Berapa banyak genom yang tidak lengkap?
2. Berapa banyak gen yang ada dalam genom?
3. Berapa lama gen yang khas?
4. Seperti apa distribusi gen di antara kromosom?

File GFF3 terbaru untuk genom manusia dapat diunduh dari sini. File README yang ada di direktori ini memberikan deskripsi singkat tentang format data ini, dan spesifikasi yang lebih lengkap dapat ditemukan di sini.

Kami akan menggunakan Pandas, komponen utama dari tumpukan SciPy yang menyediakan struktur data yang cepat, fleksibel, dan ekspresif, untuk memanipulasi dan memahami file GFF3.

Mempersiapkan

Hal pertama yang pertama, mari kita siapkan lingkungan virtual dengan tumpukan SciPy diinstal. Proses ini dapat memakan waktu jika dibangun dari sumber secara manual, karena tumpukan melibatkan banyak paket – beberapa di antaranya bergantung pada kode FORTRAN atau C eksternal. Di sini, saya sarankan menggunakan Miniconda, yang membuat pengaturannya sangat mudah.

 wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh bash Miniconda3-latest-Linux-x86_64.sh -b

Bendera -b pada baris bash memerintahkannya untuk dieksekusi dalam mode batch. Setelah perintah di atas digunakan untuk berhasil menginstal Miniconda, mulai lingkungan virtual baru untuk genomik, lalu instal tumpukan SciPy.

 mkdir -p genomics cd genomics conda create -p venv ipython matplotlib pandas

Perhatikan bahwa kami hanya menentukan 3 paket yang akan kami gunakan dalam posting ini. Jika Anda ingin semua paket terdaftar di tumpukan SciPy, cukup tambahkan ke akhir perintah conda create .

Jika Anda tidak yakin dengan nama pasti sebuah paket, coba conda search . Mari aktifkan lingkungan virtual dan mulai IPython.

 source activate venv/ ipython

IPython adalah pengganti yang jauh lebih kuat untuk antarmuka juru bahasa Python default, jadi apa pun yang biasa Anda lakukan di juru bahasa python default juga dapat dilakukan di IPython. Saya sangat merekomendasikan setiap programmer Python, yang belum pernah menggunakan IPython, untuk mencobanya.

Unduh File Anotasi

Dengan penyiapan kami sekarang selesai, mari unduh file anotasi genom manusia dalam format GFF3.

Ini adalah sekitar 37 MB, file yang sangat kecil dibandingkan dengan konten informasi genom manusia, yaitu sekitar 3 GB dalam teks biasa. Itu karena file GFF3 hanya berisi anotasi urutan, sedangkan data urutan biasanya disimpan dalam format file lain yang disebut FASTA. Jika Anda tertarik, Anda dapat mengunduh FASTA di sini, tetapi kami tidak akan menggunakan data urutan dalam tutorial ini.

 !wget ftp://ftp.ensembl.org/pub/release-85/gff3/homo_sapiens/Homo_sapiens.GRCh38.85.gff3.gz

awalan ! memberi tahu IPython bahwa ini adalah perintah Shell alih-alih perintah Python. Namun, IPython juga dapat memproses beberapa perintah shell yang sering digunakan seperti ls , pwd , rm , mkdir , rmdir bahkan tanpa awalan ! .

Perhatikan bagian atas file GFF, Anda akan melihat banyak baris metadata/pragma/direktif yang dimulai dengan ## atau #! .

Menurut README, ## berarti metadata stabil, sedangkan #! berarti itu eksperimental.

Nanti Anda juga akan melihat ### , yang merupakan arahan lain dengan makna yang lebih halus berdasarkan spesifikasi.

Komentar yang dapat dibaca manusia seharusnya setelah satu # . Untuk kesederhanaan, kami akan memperlakukan semua baris yang dimulai dengan # sebagai komentar, dan abaikan saja selama analisis kami.

 ##gff-version 3 ##sequence-region 1 1 248956422 ##sequence-region 10 1 133797422 ##sequence-region 11 1 135086622 ##sequence-region 12 1 133275309 ... ##sequence-region MT 1 16569 ##sequence-region X 1 156040895 ##sequence-region Y 2781480 56887902 #!genome-build GRCh38.p7 #!genome-version GRCh38 #!genome-date 2013-12 #!genome-build-accession NCBI:GCA_000001405.22 #!genebuild-last-updated 2016-06

Baris pertama menunjukkan bahwa versi format GFF yang digunakan dalam file ini adalah 3.

Berikut ini adalah ringkasan dari semua daerah urutan. Seperti yang akan kita lihat nanti, informasi tersebut juga dapat ditemukan di bagian tubuh file.

Garis dimulai dengan #! menunjukkan informasi tentang pembuatan genom tertentu, GRCh38.p7, yang berlaku untuk file anotasi ini.

Genome Reference Consortium (GCR) adalah konsorsium internasional, yang mengawasi pembaruan dan peningkatan beberapa rakitan genom referensi, termasuk untuk manusia, tikus, ikan zebra, dan ayam.

Memindai file ini, berikut adalah beberapa baris anotasi pertama.

 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_000001.11 ### 1 . biological_region 10469 11240 1.3e+03 . . external_name=oe %3D 0.79;logic_name=cpg 1 . biological_region 10650 10657 0.999 + . logic_name=eponine 1 . biological_region 10655 10657 0.999 - . logic_name=eponine 1 . biological_region 10678 10687 0.999 + . logic_name=eponine 1 . biological_region 10681 10688 0.999 - . logic_name=eponine ...

Kolom tersebut adalah seqid , source , type , start , end , score , strand , phase , atribut . Beberapa di antaranya sangat mudah dipahami. Ambil baris pertama sebagai contoh:

 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_000001.11

Ini adalah anotasi kromosom pertama dengan seqid 1, yang dimulai dari basis pertama hingga basis 24.895.622.

Dengan kata lain, kromosom pertama panjangnya sekitar 25 juta basa.

Analisis kami tidak memerlukan informasi dari tiga kolom dengan nilai . (yaitu skor, untai, dan fase), jadi kita bisa mengabaikannya untuk saat ini.

Kolom atribut terakhir mengatakan Chromosome 1 juga memiliki tiga nama alias, yaitu CM000663.2, chr1, dan NC_000001.11. Pada dasarnya seperti itulah tampilan file GFF3, tetapi kami tidak akan memeriksanya baris demi baris, jadi inilah saatnya memuat seluruh file ke dalam Pandas.

Pandas sangat cocok untuk menangani format GFF3 karena merupakan file tab-delimited, dan Pandas memiliki dukungan yang sangat baik untuk membaca file seperti CSV.

Perhatikan satu pengecualian untuk format tab-delimited adalah ketika GFF3 berisi ##FASTA .

Menurut spesifikasi, ##FASTA menunjukkan akhir dari bagian anotasi, yang akan diikuti dengan satu atau lebih urutan dalam format FASTA (tidak dibatasi tab). Tapi ini tidak berlaku untuk file GFF3 yang akan kita analisis.

 In [1]: import pandas as pd In [2]: pd.__version__ Out[2]: '0.18.1' In [3]: col_names = ['seqid', 'source', 'type', 'start', 'end', 'score', 'strand', 'phase', 'attributes'] Out[3]: df = pd.read_csv('Homo_sapiens.GRCh38.85.gff3.gz', compression='gzip', sep='\t', comment='#', low_memory=False, header=None, names=col_names)

Baris terakhir di atas memuat seluruh file pandas.read_csv .

Karena ini bukan file CSV standar, kami perlu sedikit menyesuaikan panggilan.

Pertama, kami memberi tahu Pandas tentang tidak tersedianya informasi header di GFF3 dengan header=None , dan kemudian kami menentukan nama yang tepat untuk setiap kolom dengan names=col_names .

Jika argumen names tidak ditentukan, Pandas akan menggunakan angka tambahan sebagai nama untuk setiap kolom.

sep='\t' memberi tahu Pandas bahwa kolom dipisahkan dengan tab, bukan dipisahkan koma. Nilai ke sep= sebenarnya bisa berupa ekspresi reguler (regex). Ini menjadi berguna jika file yang ada menggunakan pemisah yang berbeda untuk setiap kolom (oh ya, itu terjadi). comment='#' berarti baris yang dimulai dengan # dianggap sebagai komentar dan akan diabaikan.

compression='gzip' memberi tahu Pandas bahwa file input adalah file terkompresi gzip.

Selain itu, pandas.read_csv memiliki serangkaian parameter yang memungkinkan berbagai jenis format file seperti CSV untuk dibaca.

Jenis nilai yang dikembalikan adalah DataFrame , yang merupakan struktur data terpenting di Pandas, yang digunakan untuk mewakili data 2D.

Pandas juga memiliki struktur data Series dan Panel untuk data 1D dan 3D. Silakan merujuk ke dokumentasi untuk pengenalan struktur data Pandas.

Mari kita lihat beberapa entri pertama dengan metode .head .

 In [18]: df.head() Out[18]: seqid source type start end score strand phase attributes 0 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_00000... 1 1 . biological_region 10469 11240 1.3e+03 . . external_name=oe %3D 0.79;logic_name=cpg 2 1 . biological_region 10650 10657 0.999 + . logic_name=eponine 3 1 . biological_region 10655 10657 0.999 - . logic_name=eponine 4 1 . biological_region 10678 10687 0.999 + . logic_name=eponine

Outputnya diformat dengan baik dalam format tabel dengan string yang lebih panjang di kolom atribut sebagian diganti dengan ... .

Anda dapat mengatur Panda untuk tidak menghilangkan string panjang dengan pd.set_option('display.max_colwidth', -1) . Selain itu, Pandas memiliki banyak opsi yang dapat dikustomisasi.

Selanjutnya, mari kita dapatkan beberapa informasi dasar tentang kerangka data ini dengan metode .info .

 In [20]: df.info() <class 'pandas.core.frame.DataFrame'> RangeIndex: 2601849 entries, 0 to 2601848 Data columns (total 9 columns): seqid object source object type object start int64 end int64 score object strand object phase object attributes object dtypes: int64(2), object(7) memory usage: 178.7+ MB

Ini menunjukkan bahwa GFF3 memiliki 2.601.848 baris beranotasi, dan setiap baris memiliki sembilan kolom.

Untuk setiap kolom, itu juga menunjukkan tipe datanya.

Awal dan akhir itu bertipe int64, bilangan bulat yang mewakili posisi dalam genom.

Kolom lainnya semuanya bertipe object , yang mungkin berarti nilainya terdiri dari campuran bilangan bulat, float, dan string.

Ukuran semua informasi sekitar 178,7+ MB yang tersimpan di memori. Ini ternyata lebih kompak daripada file yang tidak dikompresi, yaitu sekitar 402 MB. Verifikasi cepat ditunjukkan di bawah ini.

 gunzip -c Homo_sapiens.GRCh38.85.gff3.gz > /tmp/tmp.gff3 && du -s /tmp/tmp.gff3 402M /tmp/tmp.gff3

Dari tampilan tingkat tinggi, kami telah memuat seluruh file GFF3 ke objek DataFrame dengan Python, dan semua analisis berikut kami akan didasarkan pada objek tunggal ini.

Sekarang, mari kita lihat apa itu seqid kolom pertama.

 In [29]: df.seqid.unique() Out[29]: array(['1', '10', '11', '12', '13', '14', '15', '16', '17', '18', '19', '2', '20', '21', '22', '3', '4', '5', '6', '7', '8', '9', 'GL000008.2', 'GL000009.2', 'GL000194.1', 'GL000195.1', ... 'KI270757.1', 'MT', 'X', 'Y'], dtype=object) In [30]: df.seqid.unique().shape Out[30]: (194,)

df.seqid adalah salah satu cara untuk mengakses data kolom dari kerangka data. Cara lain adalah df['seqid'] , yang merupakan sintaks yang lebih umum, karena jika nama kolom adalah kata kunci yang dicadangkan Python (misalnya class ) atau berisi . atau karakter spasi, cara pertama ( df.seqid ) tidak akan berfungsi.

Keluaran menunjukkan bahwa terdapat 194 seqid unik, yang meliputi DNA Kromosom 1 sampai 22, X, Y, dan mitokondria (MT) serta 169 seqid lainnya.

Seqid dimulai dengan KI dan GL adalah urutan DNA – atau perancah – dalam genom yang belum berhasil dirakit menjadi kromosom.

Bagi mereka yang tidak terbiasa dengan genomik, ini penting.

Meskipun rancangan genom manusia pertama keluar lebih dari 15 tahun yang lalu, genom manusia saat ini masih belum lengkap. Kesulitan dalam merakit urutan ini sebagian besar disebabkan oleh daerah berulang yang kompleks dalam genom.

Selanjutnya, mari kita lihat kolom sumber.

README mengatakan bahwa sumbernya adalah kualifikasi teks bebas yang dimaksudkan untuk menggambarkan algoritme atau prosedur operasi yang menghasilkan fitur ini.

 In [66]: df.source.value_counts() Out[66]: havana 1441093 ensembl_havana 745065 ensembl 228212 . 182510 mirbase 4701 GRCh38 194 insdc 74

Ini adalah contoh penggunaan metode value_counts , yang sangat berguna untuk penghitungan cepat variabel kategori.

Dari hasil tersebut, kita dapat melihat bahwa ada tujuh kemungkinan nilai untuk kolom ini, dan sebagian besar entri dalam file GFF3 berasal dari havana, ensembl dan ensembl_havana.

Anda dapat mempelajari lebih lanjut tentang apa arti sumber-sumber ini dan hubungan di antara mereka dalam posting ini.

Untuk mempermudah, kami akan fokus pada entri dari sumber GRCh38, havana, ensembl, dan ensembl_havan.a.

Berapa Banyak Genom yang Tidak Lengkap?

Informasi tentang setiap kromosom secara keseluruhan ada dalam entri dari sumber GRCh38, jadi mari kita menyaring sisanya, dan menetapkan hasil yang difilter ke variabel baru gdf .

 In [70]: gdf = df[df.source == 'GRCh38'] In [87]: gdf.shape Out[87]: (194, 9) In [84]: gdf.sample(10) Out[84]: seqid source type start end score strand phase attributes 2511585 KI270708.1 GRCh38 supercontig 1 127682 . . . ID=supercontig:KI270708.1;Alias=chr1_KI270708v... 2510840 GL000208.1 GRCh38 supercontig 1 92689 . . . ID=supercontig:GL000208.1;Alias=chr5_GL000208v... 990810 17 GRCh38 chromosome 1 83257441 . . . ID=chromosome:17;Alias=CM000679.2,chr17,NC_000... 2511481 KI270373.1 GRCh38 supercontig 1 1451 . . . ID=supercontig:KI270373.1;Alias=chrUn_KI270373... 2511490 KI270384.1 GRCh38 supercontig 1 1658 . . . ID=supercontig:KI270384.1;Alias=chrUn_KI270384... 2080148 6 GRCh38 chromosome 1 170805979 . . . ID=chromosome:6;Alias=CM000668.2,chr6,NC_00000... 2511504 KI270412.1 GRCh38 supercontig 1 1179 . . . ID=supercontig:KI270412.1;Alias=chrUn_KI270412... 1201561 19 GRCh38 chromosome 1 58617616 . . . ID=chromosome:19;Alias=CM000681.2,chr19,NC_000... 2511474 KI270340.1 GRCh38 supercontig 1 1428 . . . ID=supercontig:KI270340.1;Alias=chrUn_KI270340... 2594560 Y GRCh38 chromosome 2781480 56887902 . . . ID=chromosome:Y;Alias=CM000686.2,chrY,NC_00002...

Memfilter mudah di Pandas.

Jika Anda memeriksa nilai yang dievaluasi dari ekspresi df.source == 'GRCh38' , ini adalah rangkaian nilai True dan False untuk setiap entri dengan indeks yang sama dengan df . Meneruskannya ke df[] hanya akan mengembalikan entri-entri yang nilainya sesuai dengan True.

Ada 194 kunci di df[] yang df.source == 'GRCh38' .

Seperti yang telah kita lihat sebelumnya, ada juga 194 nilai unik di kolom seqid , artinya setiap entri di gdf sesuai dengan seqid tertentu.

Kemudian kami secara acak memilih 10 entri dengan metode sample untuk melihat lebih dekat.

Anda dapat melihat bahwa urutan yang belum dirakit adalah tipe supercontig sementara yang lain adalah kromosom. Untuk menghitung fraksi genom yang tidak lengkap, pertama-tama kita perlu mengetahui panjang seluruh genom, yang merupakan jumlah dari panjang semua seqid.

 In [90]: gdf = gdf.copy() In [91]: gdf['length'] = gdf.end - gdf.start + 1 In [93]: gdf.head() Out[93]: seqid source type start end score strand phase attributes length 0 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_00000... 248956421 235068 10 GRCh38 chromosome 1 133797422 . . . ID=chromosome:10;Alias=CM000672.2,chr10,NC_000... 133797421 328938 11 GRCh38 chromosome 1 135086622 . . . ID=chromosome:11;Alias=CM000673.2,chr11,NC_000... 135086621 483370 12 GRCh38 chromosome 1 133275309 . . . ID=chromosome:12;Alias=CM000674.2,chr12,NC_000... 133275308 634486 13 GRCh38 chromosome 1 114364328 . . . ID=chromosome:13;Alias=CM000675.2,chr13,NC_000... 114364327 In [97]: gdf.length.sum() Out[97]: 3096629532 In [99]: chrs = [str(_) for _ in range(1, 23)] + ['X', 'Y', 'MT'] In [101]: gdf[-gdf.seqid.isin(chrs)].length.sum() / gdf.length.sum() Out[101]: 0.0037021917421198327

Dalam cuplikan di atas terlebih dahulu, kami membuat salinan gdf dengan .copy() . Jika tidak, gdf asli hanyalah sepotong df , dan memodifikasinya secara langsung akan menghasilkan SettingWithCopyWarning (lihat di sini untuk detail lebih lanjut).

Kami kemudian menghitung panjang setiap entri dan menambahkannya kembali ke gdf sebagai kolom baru bernama "panjang". Panjang totalnya ternyata sekitar 3,1 miliar, dan fraksi dari urutan yang belum dirakit adalah sekitar 0,37%.

Berikut adalah cara kerja slicing dalam dua perintah terakhir.

Pertama, kami membuat daftar string yang mencakup semua seqid dari urutan yang dirakit dengan baik, yang semuanya adalah kromosom dan mitokondria. Kami kemudian menggunakan metode isin untuk menyaring semua entri yang seqidnya ada di daftar chrs .

Tanda minus ( - ) ditambahkan di awal indeks untuk membalikkan pemilihan, karena sebenarnya kita menginginkan semua yang tidak ada dalam daftar (yaitu kita ingin yang belum dirakit dimulai dengan KI dan GL)…

Catatan: Karena urutan yang dirakit dan tidak dirakit dibedakan berdasarkan kolom tipe, baris terakhir dapat ditulis ulang sebagai berikut untuk mendapatkan hasil yang sama.

 gdf[(gdf['type'] == 'supercontig')].length.sum() / gdf.length.sum()

Ada Berapa Gen?

Di sini kita fokus pada entri dari ensembl sumber, havana dan ensembl_havana karena merekalah tempat sebagian besar entri anotasi berada.

 In [109]: edf = df[df.source.isin(['ensembl', 'havana', 'ensembl_havana'])] In [111]: edf.sample(10) Out[111]: seqid source type start end score strand phase attributes 915996 16 havana CDS 27463541 27463592 . - 2 ID=CDS:ENSP00000457449;Parent=transcript:ENST0... 2531429 X havana exon 41196251 41196359 . + . Parent=transcript:ENST00000462850;Name=ENSE000... 1221944 19 ensembl_havana CDS 5641740 5641946 . + 0 ID=CDS:ENSP00000467423;Parent=transcript:ENST0... 243070 10 havana exon 13116267 13116340 . + . Parent=transcript:ENST00000378764;Name=ENSE000... 2413583 8 ensembl_havana exon 144359184 144359423 . + . Parent=transcript:ENST00000530047;Name=ENSE000... 2160496 6 havana exon 111322569 111322678 . - . Parent=transcript:ENST00000434009;Name=ENSE000... 839952 15 havana exon 76227713 76227897 . - . Parent=transcript:ENST00000565910;Name=ENSE000... 957782 16 ensembl_havana exon 67541653 67541782 . + . Parent=transcript:ENST00000379312;Name=ENSE000... 1632979 21 ensembl_havana exon 37840658 37840709 . - . Parent=transcript:ENST00000609713;Name=ENSE000... 1953399 4 havana exon 165464390 165464586 . + . Parent=transcript:ENST00000511992;Name=ENSE000... In [123]: edf.type.value_counts() Out[123]: exon 1180596 CDS 704604 five_prime_UTR 142387 three_prime_UTR 133938 transcript 96375 gene 42470 processed_transcript 28228 ... Name: type, dtype: int64

Metode isin digunakan lagi untuk penyaringan.

Kemudian, quick value count menunjukkan bahwa mayoritas entri adalah ekson, coding sequence (CDS), dan untranslated region (UTR).

Ini adalah elemen sub-gen, tetapi kami terutama mencari jumlah gen. Seperti yang ditunjukkan, ada 42.470, tetapi kami ingin tahu lebih banyak.

Secara khusus, siapa nama mereka, dan apa yang mereka lakukan? Untuk menjawab pertanyaan-pertanyaan tersebut, kita perlu mencermati informasi pada kolom atribut.

 In [127]: ndf = edf[edf.type == 'gene'] In [173]: ndf = ndf.copy() In [133]: ndf.sample(10).attributes.values Out[133]: array(['ID=gene:ENSG00000228611;Name=HNF4GP1;biotype=processed_pseudogene;description=hepatocyte nuclear factor 4 gamma pseudogene 1 [Source:HGNC Symbol%3BAcc:HGNC:35417];gene_id=ENSG00000228611;havana_gene=OTTHUMG00000016986;havana_version=2;logic_name=havana;version=2', 'ID=gene:ENSG00000177189;Name=RPS6KA3;biotype=protein_coding;description=ribosomal protein S6 kinase A3 [Source:HGNC Symbol%3BAcc:HGNC:10432];gene_id=ENSG00000177189;havana_gene=OTTHUMG00000021231;havana_version=5;logic_name=ensembl_havana_gene;version=12', 'ID=gene:ENSG00000231748;Name=RP11-227H15.5;biotype=antisense;gene_id=ENSG00000231748;havana_gene=OTTHUMG00000018373;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000227426;Name=VN1R33P;biotype=unitary_pseudogene;description=vomeronasal 1 receptor 33 pseudogene [Source:HGNC Symbol%3BAcc:HGNC:37353];gene_id=ENSG00000227426;havana_gene=OTTHUMG00000154474;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000087250;Name=MT3;biotype=protein_coding;description=metallothionein 3 [Source:HGNC Symbol%3BAcc:HGNC:7408];gene_id=ENSG00000087250;havana_gene=OTTHUMG00000133282;havana_version=3;logic_name=ensembl_havana_gene;version=8', 'ID=gene:ENSG00000177108;Name=ZDHHC22;biotype=protein_coding;description=zinc finger DHHC-type containing 22 [Source:HGNC Symbol%3BAcc:HGNC:20106];gene_id=ENSG00000177108;havana_gene=OTTHUMG00000171575;havana_version=3;logic_name=ensembl_havana_gene;version=5', 'ID=gene:ENSG00000249784;Name=SCARNA22;biotype=scaRNA;description=small Cajal body-specific RNA 22 [Source:HGNC Symbol%3BAcc:HGNC:32580];gene_id=ENSG00000249784;logic_name=ncrna;version=1', 'ID=gene:ENSG00000079101;Name=CLUL1;biotype=protein_coding;description=clusterin like 1 [Source:HGNC Symbol%3BAcc:HGNC:2096];gene_id=ENSG00000079101;havana_gene=OTTHUMG00000178252;havana_version=7;logic_name=ensembl_havana_gene;version=16', 'ID=gene:ENSG00000229224;Name=AC105398.3;biotype=antisense;gene_id=ENSG00000229224;havana_gene=OTTHUMG00000152025;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000255552;Name=LY6G6E;biotype=protein_coding;description=lymphocyte antigen 6 complex%2C locus G6E (pseudogene) [Source:HGNC Symbol%3BAcc:HGNC:13934];gene_id=ENSG00000255552;havana_gene=OTTHUMG00000166419;havana_version=1;logic_name=ensembl_havana_gene;version=7'], dtype=object)

Mereka diformat sebagai daftar pasangan tag-nilai yang dipisahkan titik koma. Informasi yang paling kami minati adalah nama gen, ID gen, dan deskripsi, dan kami akan mengekstraknya dengan ekspresi reguler (regex).

 import re RE_GENE_NAME = re.compile(r'Name=(?P<gene_name>.+?);') def extract_gene_name(attributes_str): res = RE_GENE_NAME.search(attributes_str) return res.group('gene_name') ndf['gene_name'] = ndf.attributes.apply(extract_gene_name)

Pertama, kami mengekstrak nama gen.

Di regex Name=(?P<gene_name>.+?); , +? digunakan sebagai ganti + karena kami ingin itu tidak serakah dan membiarkan pencarian berhenti di titik koma pertama; jika tidak, hasilnya akan sesuai dengan titik koma terakhir.

Selain itu, regex pertama kali dikompilasi dengan re.compile alih-alih digunakan secara langsung seperti dalam re.search untuk kinerja yang lebih baik karena kami akan menerapkannya ke ribuan string atribut.

extract_gene_name berfungsi sebagai fungsi pembantu untuk digunakan di pd.apply , yaitu metode yang digunakan ketika suatu fungsi perlu diterapkan pada setiap entri kerangka data atau seri.

Dalam kasus khusus ini, kami ingin mengekstrak nama gen untuk setiap entri di ndf.attributes , dan menambahkan nama kembali ke ndf di kolom baru bernama gene_name .

ID dan deskripsi gen diekstraksi dengan cara yang sama.

 RE_GENE_ID = re.compile(r'gene_id=(?P<gene_id>ENSG.+?);') def extract_gene_id(attributes_str): res = RE_GENE_ID.search(attributes_str) return res.group('gene_id') ndf['gene_id'] = ndf.attributes.apply(extract_gene_id) RE_DESC = re.compile('description=(?P<desc>.+?);') def extract_description(attributes_str): res = RE_DESC.search(attributes_str) if res is None: return '' else: return res.group('desc') ndf['desc'] = ndf.attributes.apply(extract_description)

Regex untuk RE_GENE_ID sedikit lebih spesifik karena kita tahu bahwa setiap gene_id harus dimulai dengan ENSG , di mana ENS berarti ensembl dan G berarti gen.

Untuk entri yang tidak memiliki deskripsi, kami akan mengembalikan string kosong. Setelah semuanya diekstrak, kita tidak akan menggunakan kolom atribut lagi, jadi mari kita lepaskan untuk menjaga semuanya tetap bagus dan bersih dengan metode .drop :

 In [224]: ndf.drop('attributes', axis=1, inplace=True) In [225]: ndf.head() Out[225]: seqid source type start end score strand phase gene_id gene_name desc 16 1 havana gene 11869 14409 . + . ENSG00000223972 DDX11L1 DEAD/H-box helicase 11 like 1 [Source:HGNC Sym... 28 1 havana gene 14404 29570 . - . ENSG00000227232 WASH7P WAS protein family homolog 7 pseudogene [Sourc... 71 1 havana gene 52473 53312 . + . ENSG00000268020 OR4G4P olfactory receptor family 4 subfamily G member... 74 1 havana gene 62948 63887 . + . ENSG00000240361 OR4G11P olfactory receptor family 4 subfamily G member... 77 1 ensembl_havana gene 69091 70008 . + . ENSG00000186092 OR4F5 olfactory receptor family 4 subfamily F member...

Dalam panggilan di atas, attributes menunjukkan kolom tertentu yang ingin kita jatuhkan.

axis=1 berarti kita menjatuhkan kolom alih-alih baris ( axis=0 secara default).

inplace=True berarti drop dioperasikan pada DataFrame itu sendiri alih-alih mengembalikan salinan baru dengan kolom yang ditentukan dijatuhkan.

Tampilan .head cepat menunjukkan bahwa kolom atribut memang hilang, dan tiga kolom baru: gene_name , gene_id , dan desc telah ditambahkan.

Karena penasaran, mari kita lihat apakah semua gene_id dan gene_name unik:

 In [232]: ndf.shape Out[232]: (42470, 11) In [233]: ndf.gene_id.unique().shape Out[233]: (42470,) In [234]: ndf.gene_name.unique().shape Out[234]: (42387,)

Anehnya, jumlah nama gen lebih kecil daripada ID gen, menunjukkan bahwa beberapa nama_gen harus sesuai dengan beberapa ID gen. Mari kita cari tahu apa itu.

 In [243]: count_df = ndf.groupby('gene_name').count().ix[:, 0].sort_values().ix[::-1] In [244]: count_df.head(10) Out[244]: gene_name SCARNA20 7 SCARNA16 6 SCARNA17 5 SCARNA15 4 SCARNA21 4 SCARNA11 4 Clostridiales-1 3 SCARNA4 3 C1QTNF9B-AS1 2 C11orf71 2 Name: seqid, dtype: int64 In [262]: count_df[count_df > 1].shape Out[262]: (63,) In [263]: count_df.shape Out[263]: (42387,) In [264]: count_df[count_df > 1].shape[0] / count_df.shape[0] Out[264]: 0.0014863047632528842

Kami mengelompokkan semua entri berdasarkan nilai gene_name , lalu menghitung jumlah item di setiap grup dengan .count() .

Jika Anda memeriksa keluaran dari ndf.groupby('gene_name').count() , semua kolom dihitung untuk setiap grup, tetapi sebagian besar memiliki nilai yang sama.

Perhatikan bahwa nilai NA tidak akan dipertimbangkan saat menghitung, jadi ambil saja hitungan kolom pertama, seqid ( kami menggunakan .ix[:, 0] untuk memastikan tidak ada nilai NA).

Kemudian urutkan nilai hitungan dengan .sort_values ​​dan balikkan urutannya dengan .ix[::-1] .

Hasilnya, sebuah nama gen dapat digunakan bersama hingga tujuh ID gen.

 In [255]: ndf[ndf.gene_name == 'SCARNA20'] Out[255]: seqid source type start end score strand phase gene_id gene_name desc 179399 1 ensembl gene 171768070 171768175 . + . ENSG00000253060 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 201037 1 ensembl gene 204727991 204728106 . + . ENSG00000251861 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 349203 11 ensembl gene 8555016 8555146 . + . ENSG00000252778 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 718520 14 ensembl gene 63479272 63479413 . + . ENSG00000252800 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 837233 15 ensembl gene 75121536 75121666 . - . ENSG00000252722 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 1039874 17 ensembl gene 28018770 28018907 . + . ENSG00000251818 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 1108215 17 ensembl gene 60231516 60231646 . - . ENSG00000252577 SCARNA20 small Cajal body-specific RNA 20 [Source:HGNC Symbol%3BAcc:HGNC:32578]

Melihat lebih dekat pada semua gen SCARNA20 menunjukkan bahwa mereka memang berbeda.

Meskipun mereka memiliki nama yang sama, mereka berada di posisi genom yang berbeda.

Namun, deskripsi mereka tampaknya tidak terlalu membantu dalam membedakan mereka.

Intinya di sini adalah untuk mengetahui bahwa nama gen tidak unik untuk semua ID gen, dan sekitar 0,15% dari nama itu dimiliki oleh banyak gen.

Berapa Lama Gen Khas?

Mirip dengan apa yang kami lakukan ketika kami mencoba memahami ketidaklengkapan genom, kami dapat dengan mudah menambahkan kolom length ke ndf :

 In [277]: ndf['length'] = ndf.end - ndf.start + 1 In [278]: ndf.length.describe() Out[278]: count 4.247000e+04 mean 3.583348e+04 std 9.683485e+04 min 8.000000e+00 25% 8.840000e+02 50% 5.170500e+03 75% 3.055200e+04 max 2.304997e+06 Name: length, dtype: float64

.describe() menghitung beberapa statistik sederhana berdasarkan nilai panjang:

• Panjang rata-rata gen adalah sekitar 36.000 basa

• Panjang rata-rata gen adalah sekitar 5.200 basa

• Panjang gen minimum dan maksimum masing-masing sekitar delapan dan 2,3 juta basa.

Karena mean jauh lebih besar dari median, ini menyiratkan bahwa distribusi panjang miring ke kanan. Untuk melihat lebih konkret, mari kita plot distribusinya.

 import matplotlib as plt ndf.length.plot(kind='hist', bins=50, logy=True) plt.show()

Pandas menyediakan antarmuka sederhana untuk matplotlib untuk membuat plot sangat berguna dengan DataFrames atau seri.

Dalam hal ini, dikatakan bahwa kita menginginkan plot histogram ( kind='hist' ) dengan 50 bin, dan biarkan sumbu y berada pada skala log ( logy=True ).

Dari histogram, kita dapat melihat bahwa sebagian besar gen berada di dalam bin pertama. Namun, beberapa panjang gen bisa lebih dari dua juta basa. Mari kita cari tahu apa itu:

 In [39]: ndf[ndf.length > 2e6].sort_values('length').ix[::-1] Out[39]: seqid source type start end score strand phase gene_name gene_id desc length 2309345 7 ensembl_havana gene 146116002 148420998 . + . CNTNAP2 ENSG00000174469 contactin associated protein-like 2 [Source:HG... 2304997 2422510 9 ensembl_havana gene 8314246 10612723 . - . PTPRD ENSG00000153707 protein tyrosine phosphatase%2C receptor type ... 2298478 2527169 X ensembl_havana gene 31097677 33339441 . - . DMD ENSG00000198947 dystrophin [Source:HGNC Symbol%3BAcc:HGNC:2928] 2241765 440886 11 ensembl_havana gene 83455012 85627922 . - . DLG2 ENSG00000150672 discs large MAGUK scaffold protein 2 [Source:H... 2172911 2323457 8 ensembl_havana gene 2935353 4994972 . - . CSMD1 ENSG00000183117 CUB and Sushi multiple domains 1 [Source:HGNC ... 2059620 1569914 20 ensembl_havana gene 13995369 16053197 . + . MACROD2 ENSG00000172264 MACRO domain containing 2 [Source:HGNC Symbol%... 2057829

Seperti yang Anda lihat, gen terpanjang bernama CNTNAP2, yang merupakan kependekan dari contactin linked protein-like 2. Menurut halaman wikipedia-nya,

Gen ini mencakup hampir 1,6% dari kromosom 7 dan merupakan salah satu gen terbesar dalam genom manusia.

Memang! Kami baru saja memverifikasinya sendiri. Sebaliknya, bagaimana dengan gen terkecil? Ternyata mereka bisa sesingkat delapan pangkalan.

 In [40]: ndf.sort_values('length').head() Out[40]: seqid source type start end score strand phase gene_name gene_id desc length 682278 14 havana gene 22438547 22438554 . + . TRDD1 ENSG00000223997 T cell receptor delta diversity 1 [Source:HGNC... 8 682282 14 havana gene 22439007 22439015 . + . TRDD2 ENSG00000237235 T cell receptor delta diversity 2 [Source:HGNC... 9 2306836 7 havana gene 142786213 142786224 . + . TRBD1 ENSG00000282431 T cell receptor beta diversity 1 [Source:HGNC ... 12 682286 14 havana gene 22449113 22449125 . + . TRDD3 ENSG00000228985 T cell receptor delta diversity 3 [Source:HGNC... 13 1879625 4 havana gene 10238213 10238235 . - . AC006499.9 ENSG00000271544 23

Panjang dari dua kasus ekstrim tersebut adalah lima orde besarnya terpisah (2,3 juta vs 8), yang sangat besar dan yang dapat menjadi indikasi tingkat keragaman kehidupan.

Sebuah gen tunggal dapat diterjemahkan ke banyak protein berbeda melalui proses yang disebut penyambungan alternatif, sesuatu yang belum kami jelajahi. Informasi tersebut juga ada di dalam file GFF3, tetapi di luar cakupan posting ini.

Distribusi Gen Antara Kromosom

Hal terakhir yang ingin saya diskusikan adalah distribusi gen di antara kromosom, yang juga berfungsi sebagai contoh untuk memperkenalkan metode .merge untuk menggabungkan dua DataFrames. Secara intuitif, kromosom yang lebih panjang kemungkinan menampung lebih banyak gen. Mari kita lihat apakah itu benar.

 In [53]: ndf = ndf[ndf.seqid.isin(chrs)] In [54]: chr_gene_counts = ndf.groupby('seqid').count().ix[:, 0].sort_values().ix[::-1] Out[54]: chr_gene_counts seqid 1 3902 2 2806 11 2561 19 2412 17 2280 3 2204 6 2154 12 2140 7 2106 5 2002 16 1881 X 1852 4 1751 9 1659 8 1628 10 1600 15 1476 14 1449 22 996 20 965 13 872 18 766 21 541 Y 436 Name: source, dtype: int64

Kami meminjam variabel chrs dari bagian sebelumnya, dan menggunakannya untuk menyaring urutan yang belum dirakit. Berdasarkan keluarannya, Kromosom 1 terbesar memang memiliki gen paling banyak. Sementara Kromosom Y memiliki jumlah gen terkecil, itu bukan kromosom terkecil.

Perhatikan bahwa tampaknya tidak ada gen dalam mitokondria (MT), yang tidak benar.

Sedikit lebih banyak penyaringan pada DataFrame df pertama yang dikembalikan oleh pd.read_csv menunjukkan bahwa semua gen MT berasal dari sumber insdc (yang disaring sebelumnya saat membuat edf di mana kami hanya mempertimbangkan sumber havana, ensembl, atau ensembl_havana).

 In [134]: df[(df.type == 'gene') & (df.seqid == 'MT')] Out[134]: seqid source type start end score strand phase attributes 2514003 MT insdc gene 648 1601 . + . ID=gene:ENSG00000211459;Name=MT-RNR1;biotype=M... 2514009 MT insdc gene 1671 3229 . + . ID=gene:ENSG00000210082;Name=MT-RNR2;biotype=M... 2514016 MT insdc gene 3307 4262 . + . ID=gene:ENSG00000198888;Name=MT-ND1;biotype=pr... 2514029 MT insdc gene 4470 5511 . + . ID=gene:ENSG00000198763;Name=MT-ND2;biotype=pr... 2514048 MT insdc gene 5904 7445 . + . ID=gene:ENSG00000198804;Name=MT-CO1;biotype=pr... 2514058 MT insdc gene 7586 8269 . + . ID=gene:ENSG00000198712;Name=MT-CO2;biotype=pr... 2514065 MT insdc gene 8366 8572 . + . ID=gene:ENSG00000228253;Name=MT-ATP8;biotype=p... 2514069 MT insdc gene 8527 9207 . + . ID=gene:ENSG00000198899;Name=MT-ATP6;biotype=p... 2514073 MT insdc gene 9207 9990 . + . ID=gene:ENSG00000198938;Name=MT-CO3;biotype=pr... 2514080 MT insdc gene 10059 10404 . + . ID=gene:ENSG00000198840;Name=MT-ND3;biotype=pr... 2514087 MT insdc gene 10470 10766 . + . ID=gene:ENSG00000212907;Name=MT-ND4L;biotype=p... 2514091 MT insdc gene 10760 12137 . + . ID=gene:ENSG00000198886;Name=MT-ND4;biotype=pr... 2514104 MT insdc gene 12337 14148 . + . ID=gene:ENSG00000198786;Name=MT-ND5;biotype=pr... 2514108 MT insdc gene 14149 14673 . - . ID=gene:ENSG00000198695;Name=MT-ND6;biotype=pr... 2514115 MT insdc gene 14747 15887 . + . ID=gene:ENSG00000198727;Name=MT-CYB;biotype=pr...

This example also shows how to combine two conditions during filtering with & ; the logical operator for “or” would be | .

Note that the parentheses around each condition are required, and this part of the syntax in Pandas is different from Python, which would have been literal and and or .

Next, let's borrow the gdf DataFrame from the previous section as a source for the length of each chromosome:

 In [61]: gdf = gdf[gdf.seqid.isin(chrs)] In [62]: gdf.drop(['start', 'end', 'score', 'strand', 'phase' ,'attributes'], axis=1, inplace=True) In [63]: gdf.sort_values('length').ix[::-1] Out[63]: seqid source type length 0 1 GRCh38 chromosome 248956422 1364641 2 GRCh38 chromosome 242193529 1705855 3 GRCh38 chromosome 198295559 1864567 4 GRCh38 chromosome 190214555 1964921 5 GRCh38 chromosome 181538259 2080148 6 GRCh38 chromosome 170805979 2196981 7 GRCh38 chromosome 159345973 2514125 X GRCh38 chromosome 156040895 2321361 8 GRCh38 chromosome 145138636 2416560 9 GRCh38 chromosome 138394717 328938 11 GRCh38 chromosome 135086622 235068 10 GRCh38 chromosome 133797422 483370 12 GRCh38 chromosome 133275309 634486 13 GRCh38 chromosome 114364328 674767 14 GRCh38 chromosome 107043718 767312 15 GRCh38 chromosome 101991189 865053 16 GRCh38 chromosome 90338345 990810 17 GRCh38 chromosome 83257441 1155977 18 GRCh38 chromosome 80373285 1559144 20 GRCh38 chromosome 64444167 1201561 19 GRCh38 chromosome 58617616 2594560 Y GRCh38 chromosome 54106423 1647482 22 GRCh38 chromosome 50818468 1616710 21 GRCh38 chromosome 46709983 2513999 MT GRCh38 chromosome 16569

The columns that are not relevant to the analysis are dropped for clarity.

Yes, .drop can also take a list of columns and drop them altogether in one operation.

Note that the row with a seqid of MT is still there; we will get back to it later. The next operation we will perform is merge the two datasets based on the values of seqid.

 In [73]: cdf = chr_gene_counts.to_frame(name='gene_count').reset_index() In [75]: cdf.head(2) Out[75]: seqid gene_count 0 1 3902 1 2 2806 In [78]: merged = gdf.merge(cdf, on='seqid') In [79]: merged Out[79]: seqid source type length gene_count 0 1 GRCh38 chromosome 248956422 3902 1 10 GRCh38 chromosome 133797422 1600 2 11 GRCh38 chromosome 135086622 2561 3 12 GRCh38 chromosome 133275309 2140 4 13 GRCh38 chromosome 114364328 872 5 14 GRCh38 chromosome 107043718 1449 6 15 GRCh38 chromosome 101991189 1476 7 16 GRCh38 chromosome 90338345 1881 8 17 GRCh38 chromosome 83257441 2280 9 18 GRCh38 chromosome 80373285 766 10 19 GRCh38 chromosome 58617616 2412 11 2 GRCh38 chromosome 242193529 2806 12 20 GRCh38 chromosome 64444167 965 13 21 GRCh38 chromosome 46709983 541 14 22 GRCh38 chromosome 50818468 996 15 3 GRCh38 chromosome 198295559 2204 16 4 GRCh38 chromosome 190214555 1751 17 5 GRCh38 chromosome 181538259 2002 18 6 GRCh38 chromosome 170805979 2154 19 7 GRCh38 chromosome 159345973 2106 20 8 GRCh38 chromosome 145138636 1628 21 9 GRCh38 chromosome 138394717 1659 22 X GRCh38 chromosome 156040895 1852 23 Y GRCh38 chromosome 54106423 436

Since chr_gene_counts is still a Series object, which doesn't support a merge operation, we need to convert it to a DataFrame object first with .to_frame .

.reset_index() converts the original index (ie seqid ) into a new column and resets current index as 0-based incremental numbers.

The output from cdf.head(2) shows what it looks like. Next, we used the .merge method to combine the two DataFrame on the seqid column ( on='seqid' ).

After merging gdf and cdf , the MT entry is still missing. This is because, by default, .merge operates an inner join, while left join, right join, or outer join are available by tuning the how parameter.

Please refer to the documentation for more details.

Later, you may find that there is also a related .join method. .merge and .join are similar yet have different APIs.

According to the official documentation says

The related DataFrame.join method, uses merge internally for the index-on-index and index-on-column(s) joins, but joins on indexes by default rather than trying to join on common columns (the default behavior for merge). If you are joining on index, you may wish to use DataFrame.join to save yourself some typing.

Basically, .merge is more general-purpose and used by .join .

Finally, we are ready to calculate the correlation between chromosome length and gene_count .

 In [81]: merged[['length', 'gene_count']].corr() Out[81]: length gene_count length 1.000000 0.728221 gene_count 0.728221 1.000000

By default .corr calculates the Pearson correlation between all pairs of columns in a dataframe.

But we have only a single pair of columns in this case, and the correlation turns out to be positive – 0.73.

In other words, the larger the chromosome, the more likely it is to have more genes. Let's also plot the two columns after sorting the value pairs by length :

 ax = merged[['length', 'gene_count']].sort_values('length').plot(x='length', y='gene_count',) # add some margin to both ends of x axis xlim = ax.get_xlim() margin = xlim[0] * 0.1 ax.set_xlim([xlim[0] - margin, xlim[1] + margin]) # Label each point on the graph for (s, x, y) in merged[['seqid', 'length', 'gene_count']].sort_values('length').values: ax.text(x, y - 100, str(s)) 

As seen in image above, even though it is a positive correlation overall, it does not hold for all chromosomes. In particular, for Chromosome 17, 16, 15, 14, 13, the correlation is actually negative, meaning the number of genes on the chromosome decreases as the chromosome size increases.

Findings and Future Research

That ends our tutorial on the manipulation of an annotation file for human genome in GFF3 format with the SciPy stack. The tools we've mainly used include IPython, Pandas, and matplotlib. During the tutorial, not only have we learned some of the most common and useful operations in Pandas, we also answered some very interesting questions about our genome. In summary:

1. Sekitar 0,37% dari genom manusia masih belum lengkap meskipun draf pertama keluar lebih dari 15 tahun yang lalu.
2. Ada sekitar 42.000 gen dalam genom manusia berdasarkan file GFF3 khusus yang kami gunakan ini.
3. Panjang gen dapat berkisar dari beberapa lusin hingga lebih dari dua juta basa.
4. Gen tidak terdistribusi secara merata di antara kromosom. Secara keseluruhan, semakin besar kromosom, semakin banyak gen yang ditampungnya, tetapi untuk subset kromosom, korelasinya bisa negatif.

File GFF3 sangat kaya akan informasi anotasi, dan kami baru saja menggores permukaannya. Jika Anda tertarik untuk eksplorasi lebih lanjut, berikut adalah beberapa pertanyaan yang dapat Anda mainkan:

1. Berapa banyak transkrip yang biasanya dimiliki gen? Berapa persentase gen yang memiliki lebih dari 1 transkrip?
2. Berapa banyak isoform yang biasanya dimiliki transkrip?
3. Berapa banyak ekson, CDS, dan UTR yang biasanya dimiliki transkrip? Apa ukuran mereka?
4. Apakah mungkin untuk mengkategorikan gen berdasarkan fungsinya seperti yang dijelaskan pada kolom deskripsi?