SciPyスタックを使用したゲノムの包括的な紹介

公開: 2022-03-11

バイオインフォマティクスは、生物学的データを分析および理解するための方法とソフトウェアツールを開発する学際的な分野です。

もっと簡単に言えば、生物学のデータサイエンスと簡単に考えることができます。

多くの種類の生物学的データの中で、ゲノミクスデータは最も広く分析されているものの1つです。特に次世代DNAシーケンシング（NGS）テクノロジーの急速な進歩に伴い、ゲノミクスデータの量は飛躍的に増加しています。 Stephens、Zachary D et al。によると、ゲノミクスデータの取得は1年あたりのエクサバイト規模です。

SciPyは、科学計算用の多くのPythonモジュールを収集します。これは、多くのバイオインフォマティクスのニーズに理想的です。

つぶやき

この投稿では、SciPyスタックを使用してヒトゲノムのGFF3ファイルを分析する例をデモします。 Generic Feature Formatバージョン3（GFF3）は、ゲノム機能を保存するための現在の標準テキストファイル形式です。特に、この投稿では、SciPyスタックを使用してヒトゲノムに関する次の質問に答える方法を学習します。

ゲノムのどのくらいが不完全ですか？
ゲノムにはいくつの遺伝子がありますか？
典型的な遺伝子の長さはどれくらいですか？
染色体間の遺伝子分布はどのように見えますか？

ヒトゲノムの最新のGFF3ファイルはこちらからダウンロードできます。このディレクトリにあるREADMEファイルには、このデータ形式の簡単な説明が記載されており、より詳細な仕様がここにあります。

SciPyスタックの主要コンポーネントであるPandasを使用して、高速で柔軟で表現力豊かなデータ構造を提供し、GFF3ファイルを操作および理解します。

設定

まず最初に、SciPyスタックがインストールされた仮想環境をセットアップしましょう。スタックには多くのパッケージが含まれるため、ソースから手動でビルドする場合、このプロセスには時間がかかる可能性があります。その一部は外部のFORTRANまたはCコードに依存します。ここでは、セットアップが非常に簡単なMinicondaの使用をお勧めします。

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

bash行の-bフラグは、バッチモードで実行するように指示します。上記のコマンドを使用してMinicondaを正常にインストールした後、ゲノミクス用の新しい仮想環境を開始してから、SciPyスタックをインストールします。

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

この投稿で使用するパッケージは3つだけ指定していることに注意してください。 SciPyスタックにリストされているすべてのパッケージが必要な場合は、それらをconda createコマンドの最後に追加するだけです。

パッケージの正確な名前がわからない場合は、 conda searchを試してください。仮想環境をアクティブにして、IPythonを起動しましょう。

 source activate venv/ ipython

IPythonは、デフォルトのPythonインタープリターインターフェイスに代わる非常に強力なものであるため、デフォルトのPythonインタープリターで行っていた処理はすべてIPythonでも実行できます。まだIPythonを使用していないすべてのPythonプログラマーに、試してみることを強くお勧めします。

注釈ファイルをダウンロードする

セットアップが完了したら、GFF3形式のヒトゲノムアノテーションファイルをダウンロードしましょう。

これは約37MBであり、プレーンテキストで約3GBであるヒトゲノムの情報コンテンツと比較して非常に小さいファイルです。これは、GFF3ファイルにはシーケンスの注釈のみが含まれているのに対し、シーケンスデータは通常FASTAと呼ばれる別のファイル形式で保存されているためです。興味のある方は、ここからFASTAをダウンロードできますが、このチュートリアルではシーケンスデータを使用しません。

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

接頭辞! これがPythonコマンドではなくシェルコマンドであることをIPythonに通知します。ただし、IPythonは、接頭辞が付いていなくても、 ls 、 pwd 、 rm 、 mkdir 、 rmdirなどの頻繁に使用されるシェルコマンドを処理することもできます! 。

GFFファイルの先頭を見ると、 ##または#!で始まる多くのメタデータ/プラグマ/ディレクティブの行が表示されます。。

READMEによると、 ##はメタデータが安定していることを意味し、 #! それは実験的であることを意味します。

後で、 ###も表示されます。これは、仕様に基づいてさらに微妙な意味を持つ別のディレクティブです。

人間が読めるコメントは、単一の#の後にあるはずです。簡単にするために、 #で始まるすべての行をコメントとして扱い、分析中は無視します。

 ##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

最初の行は、このファイルで使用されているGFF形式のバージョンが3であることを示しています。

続いて、すべてのシーケンス領域の要約があります。後で説明するように、このような情報はファイルの本文部分にもあります。

#! は、このアノテーションファイルが適用されるゲノムの特定のビルドGRCh38.p7に関する情報を示しています。

ゲノムリファレンスコンソーシアム（GCR）は国際的なコンソーシアムであり、ヒト、マウス、ゼブラフィッシュ、ニワトリなど、いくつかのリファレンスゲノムアセンブリの更新と改善を監督しています。

このファイルをスキャンすると、最初の数行の注釈が表示されます。

 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 ...

列は、 seqid 、 source 、 type 、 start 、 end 、 score 、 strand 、 phase 、 attributesです。それらのいくつかは非常に理解しやすいです。例として最初の行を取り上げます。

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

これは、1番染色体から24,895,622番目の塩基までの1番染色体の注釈です。

言い換えれば、最初の染色体は約2500万塩基の長さです。

分析では、値が。の3つの列からの情報は必要ありません. （つまり、スコア、ストランド、フェーズ）なので、今のところは無視できます。

最後の属性列には、1番染色体にも3つのエイリアス名、つまりCM000663.2、chr1、およびNC_000001.11があることが示されています。これは基本的にGFF3ファイルの外観ですが、行ごとに検査することはないので、ファイル全体をPandasにロードするときが来ました。

Pandasはタブ区切りファイルであるためGFF3形式の処理に最適であり、PandasはCSVのようなファイルの読み取りを非常によくサポートしています。

タブ区切り形式の1つの例外は、GFF3に##FASTAが含まれている場合です。

仕様によれば、 ##FASTAは注釈部分の終わりを示し、その後にFASTA（タブ区切りではない）形式の1つ以上のシーケンスが続きます。ただし、これは、分析するGFF3ファイルには当てはまりません。

 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)

上記の最後の行は、 pandas.read_csvメソッドを使用してGFF3ファイル全体をロードします。

これは標準のCSVファイルではないため、呼び出しを少しカスタマイズする必要があります。

まず、 header=NoneでGFF3のヘッダー情報が利用できないことをパンダに通知し、次にnames=col_namesで各列の正確な名前を指定します。

names引数が指定されていない場合、Pandasは各列の名前として増分番号を使用します。

sep='\t'は、列がコンマ区切りではなくタブ区切りであることをパンダに通知します。 sep=の値は、実際には正規表現（regex）にすることができます。これは、手元のファイルが列ごとに異なる区切り文字を使用している場合に便利です（そうそう、そうなります）。 comment='#'は、 #で始まる行がコメントと見なされ、無視されることを意味します。

compression='gzip'は、入力ファイルがgzip圧縮されていることをPandasに通知します。

さらに、 pandas.read_csvには、さまざまな種類のCSVのようなファイル形式を読み取ることができる豊富なパラメータセットがあります。

戻り値のタイプはDataFrame 。これは、2Dデータを表すために使用されるPandasで最も重要なデータ構造です。

Pandasには、それぞれ1Dおよび3Dデータ用のSeriesおよびPanelデータ構造もあります。 Pandasのデータ構造の概要については、ドキュメントを参照してください。

.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

出力は表形式で適切にフォーマットされ、属性列の長い文字列が部分的に...に置き換えられています。

pd.set_option('display.max_colwidth', -1)を使用して、長い文字列を省略しないようにPandasを設定できます。さらに、パンダにはカスタマイズ可能な多くのオプションがあります。

次に、 .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

これは、GFF3に2,601,848の注釈付き行があり、各行に9つの列があることを示しています。

各列について、データ型も表示されます。

その開始と終了はint64タイプであり、整数はゲノム内の位置を表します。

他の列はすべてobject型です。これは、おそらくそれらの値が整数、浮動小数点数、および文字列の混合で構成されていることを意味します。

すべての情報のサイズは、メモリに保存されている約178.7MB以上です。これは、約402MBの非圧縮ファイルよりもコンパクトであることがわかります。簡単な検証を以下に示します。

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

高レベルの観点から、GFF3ファイル全体をPythonのDataFrameオブジェクトにロードしました。以降の分析はすべて、この単一のオブジェクトに基づいています。

それでは、最初の列seqidが何であるかを見てみましょう。

 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は、データフレームから列データにアクセスする1つの方法です。もう1つの方法は、 df['seqid']です。これは、より一般的な構文です。これは、列名がPythonの予約キーワード（例： class ）であるか、。が含まれているため. またはスペース文字の場合、最初の方法（ df.seqid ）は機能しません。

出力は、染色体1から22、X、Y、およびミトコンドリア（MT）DNAを含む194の固有の配列と、169の他の配列があることを示しています。

KIおよびGLで始まるseqidは、染色体に正常に組み立てられていないゲノム内のDNA配列（または足場）です。

ゲノミクスに不慣れな人にとって、これは重要です。

最初のヒトゲノムドラフトは15年以上前に発表されましたが、現在のヒトゲノムはまだ不完全です。これらの配列を組み立てるのが難しいのは、主にゲノム内の複雑な反復領域が原因です。

次に、ソース列を見てみましょう。

READMEによると、ソースは、この機能を生成したアルゴリズムまたは操作手順を説明することを目的としたフリーテキスト修飾子です。

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

これは、 value_countsメソッドの使用例です。これは、カテゴリ変数のクイックカウントに非常に役立ちます。

結果から、この列には7つの可能な値があり、GFF3ファイルのエントリの大部分はhavana、ensembl、およびensembl_havanaからのものであることがわかります。

この投稿で、これらのソースの意味とそれらの間の関係について詳しく知ることができます。

物事を単純にするために、ソースGRCh38、havana、ensembl、およびensembl_havan.aからのエントリに焦点を当てます。

ゲノムのどのくらいが不完全ですか？

各染色体全体に関する情報はソースGRCh38からのエントリにあるので、最初に残りを除外し、フィルタリングされた結果を新しい変数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...

パンダではフィルタリングが簡単です。

式df.source == 'GRCh38'から評価された値を調べると、これはdfと同じインデックスを持つ各エントリの一連のTrue値とFalse値です。これをdf[]に渡すと、対応する値がTrueであるエントリのみが返されます。

df.source == 'GRCh38'であるdf[]には194個のキーがあります。

前に見たように、 seqid列にも194の一意の値があります。これは、gdfの各エントリが特定のgdfに対応することを意味します。

次に、 sampleメソッドを使用してランダムに10個のエントリを選択し、詳しく調べます。

組み立てられていないシーケンスはスーパーコンティグタイプであり、他のシーケンスは染色体であることがわかります。不完全なゲノムの割合を計算するには、最初にゲノム全体の長さを知る必要があります。これは、すべての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

上記のスニペットでは、最初に.copy()を使用してgdfのコピーを作成しました。それ以外の場合、元のgdfはdfのスライスにすぎず、直接変更すると、 SettingWithCopyWarningになります（詳細については、ここを参照してください）。

次に、各エントリの長さを計算し、「length」という名前の新しい列としてgdfに追加し直します。全長は約31億であり、アセンブルされていない配列の割合は約0.37％です。

最後の2つのコマンドでスライスがどのように機能するかを次に示します。

まず、すべての染色体とミトコンドリアである、よく組み立てられたシーケンスのすべてのシーケンスをカバーする文字列のリストを作成します。次に、 isinメソッドを使用して、seqidがchrsリストにあるすべてのエントリをフィルタリングします。

インデックスの先頭にマイナス記号（ - ）を追加して、選択を逆にします。これは、リストにないものすべてが実際に必要なためです（つまり、KIとGLで始まるアセンブルされていないものが必要です）…

注：アセンブルされたシーケンスとアセンブルされていないシーケンスはタイプ列によって区別されるため、最後の行を次のように書き換えて同じ結果を得ることができます。

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

遺伝子はいくつありますか？

ここでは、ソースensembl、havana、およびensembl_havanaからのエントリに焦点を当てます。これらは、アノテーションエントリの大部分が属する場所だからです。

 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

isinメソッドはフィルタリングに再び使用されます。

次に、値をすばやくカウントすると、エントリの大部分がエクソン、コーディングシーケンス（CDS）、および非翻訳領域（UTR）であることがわかります。

これらはサブ遺伝子要素ですが、私たちは主に遺伝子数を探しています。示されているように、42,470がありますが、もっと知りたいです。

具体的には、彼らの名前は何ですか、そして彼らは何をしますか？これらの質問に答えるには、属性列の情報を詳しく調べる必要があります。

 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)

これらは、タグと値のペアのセミコロンで区切られたリストとしてフォーマットされます。私たちが最も興味を持っている情報は、遺伝子名、遺伝子ID、説明であり、正規表現（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)

まず、遺伝子名を抽出します。

正規表現でName=(?P<gene_name>.+?); 、 +? +の代わりに使用されます。これは、貪欲ではなく、最初のセミコロンで検索を停止するためです。それ以外の場合、結果は最後のセミコロンまで一致します。

また、正規表現は、数千の属性文字列に適用されるため、パフォーマンスを向上させるためにre.searchのように直接使用されるのではなく、最初にre.compileでコンパイルされます。

extract_gene_nameは、 pd.applyで使用されるヘルパー関数として機能します。これは、データフレームまたはシリーズのすべてのエントリに関数を適用する必要がある場合に使用するメソッドです。

この特定のケースでは、 ndf.attributesのすべてのエントリの遺伝子名を抽出し、 ndfという新しい列のgene_nameに名前を追加し直します。

遺伝子IDと説明は同様の方法で抽出されます。

 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)

RE_GENE_IDの正規表現は、すべてのgene_idがENSGで始まる必要があることがわかっているため、もう少し具体的です。ここで、 ENSはensemblを意味し、 Gは遺伝子を意味します。

説明がないエントリの場合は、空の文字列を返します。すべてが抽出された後、属性列は使用しなくなります。そのため、メソッド.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...

上記の呼び出しでは、 attributesは削除する特定の列を示しています。

axis=1は、行ではなく列を削除することを意味します（デフォルトではaxis=0 ）。

inplace=Trueは、指定された列がドロップされた新しいコピーを返すのではなく、DataFrame自体でドロップが操作されることを意味します。

.headをすばやく確認すると、属性列が実際になくなっており、 gene_name 、 gene_id 、 descの3つの新しい列が追加されていることがわかります。

好奇心から、すべてのgene_idとgene_nameが一意であるかどうかを見てみましょう。

 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,)

驚いたことに、遺伝子名の数は遺伝子IDの数よりも少なく、いくつかのgene_nameが複数の遺伝子IDに対応している必要があることを示しています。それらが何であるかを調べましょう。

 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

すべてのエントリをgene_nameの値でグループ化し、 .count()を使用して各グループのアイテム数をカウントします。

ndf.groupby('gene_name').count()からの出力を調べると、すべての列がグループごとにカウントされますが、ほとんどの列の値は同じです。

NA値はカウント時に考慮されないため、最初の列seqidのカウントのみを取得することに注意してください（NA値がないことを確認するために.ix[:, 0]を使用します）。

次に、カウント値を.sort_valuesで並べ替え、 .ix[::-1]で順序を逆にします。

その結果、遺伝子名は最大7つの遺伝子IDと共有できます。

 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]

すべてのSCARNA20遺伝子を詳しく見ると、実際にはすべて異なっていることがわかります。

それらは同じ名前を共有していますが、ゲノムの異なる位置にあります。

ただし、それらの説明は、それらを区別するのにあまり役立たないようです。

ここでのポイントは、遺伝子名がすべての遺伝子IDで一意であるとは限らず、複数の遺伝子で共有されている名前の約0.15％であることを知っておく必要があります。

典型的な遺伝子の長さはどれくらいですか？

ゲノムの不完全性を理解しようとしたときに行ったことと同様に、 ndfにlengthの列を簡単に追加できます。

 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()は、長さの値に基づいていくつかの簡単な統計を計算します。

遺伝子の平均の長さは約36,000塩基です
遺伝子の長さの中央値は約5,200塩基長です
最小遺伝子長と最大遺伝子長は、それぞれ約800万塩基と230万塩基の長さです。

平均は中央値よりもはるかに大きいため、長さの分布が右に偏っていることを意味します。より具体的な外観にするために、分布をプロットしてみましょう。

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

Pandasは、matplotlibへのシンプルなインターフェイスを提供し、DataFrameまたはシリーズでのプロットを非常に便利にします。

この場合、50個のビンを持つヒストグラムプロット（ kind='hist' ）が必要であり、y軸を対数スケール（ logy=True ）にすることを示しています。

ヒストグラムから、遺伝子の大部分が最初のビン内にあることがわかります。ただし、一部の遺伝子の長さは200万塩基を超える場合があります。それらが何であるかを調べましょう：

 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

ご覧のとおり、最長の遺伝子はCNTNAP2という名前で、これはコンタクチン関連タンパク質様2の略です。そのウィキペディアのページによると、

この遺伝子は7番染色体のほぼ1.6％を含み、ヒトゲノムで最大の遺伝子の1つです。

それはそう！私たちは自分たち自身を検証しました。対照的に、最小の遺伝子はどうですか？それらは8塩基まで短くなる可能性があることがわかりました。

 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

2つの極端なケースの長さは5桁離れており（230万対8）、これは非常に大きく、生命の多様性のレベルを示している可能性があります。

単一の遺伝子は、選択的スプライシングと呼ばれるプロセスを介して多くの異なるタンパク質に翻訳することができますが、これは私たちが調査していないことです。このような情報もGFF3ファイル内にありますが、この投稿の範囲外です。

染色体間の遺伝子分布

最後に説明したいのは、染色体間の遺伝子分布です。これは、2つのDataFrameを組み合わせるための.mergeメソッドを導入するための例としても役立ちます。直感的には、染色体が長いほど、より多くの遺伝子をホストする可能性があります。それが本当かどうか見てみましょう。

 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

前のセクションからchrs変数を借用し、それを使用して、アセンブルされていないシーケンスを除外しました。出力に基づいて、最大の染色体1は確かに最も多くの遺伝子を持っています。染色体Yは遺伝子の数が最も少ないですが、それは最小の染色体ではありません。

ミトコンドリア（MT）には遺伝子がないように見えることに注意してください。これは真実ではありません。

dfによって返される最初のpd.read_csvをもう少しフィルタリングすると、すべてのMT遺伝子がソースinsdc（havana、ensembl、またはensembl_havanaのソースのみを考慮したedfを生成する前にフィルターで除外された）からのものであることがわかります。

 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. 要約すれば：

15年以上前に最初のドラフトが発表されたにもかかわらず、ヒトゲノムの約0.37％はまだ不完全です。
私たちが使用したこの特定のGFF3ファイルに基づくと、ヒトゲノムには約42,000の遺伝子があります。
遺伝子の長さは、数十から200万塩基以上の範囲です。
遺伝子は染色体間で均等に分布していません。全体として、染色体が大きいほど、それがホストする遺伝子は多くなりますが、染色体のサブセットの場合、相関は負になる可能性があります。

GFF3ファイルは注釈情報が非常に豊富で、表面をかじったところです。さらに詳しく調べることに興味がある場合は、次の質問を試してみてください。

遺伝子には通常いくつの転写産物がありますか？遺伝子の何パーセントが複数の転写産物を持っていますか？
トランスクリプトには通常、いくつのアイソフォームがありますか？
転写産物には通常、いくつのエクソン、CDS、およびUTRがありますか？それらはどのようなサイズですか？
説明欄に記載されているように、機能に基づいて遺伝子を分類することは可能ですか？