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

もっと簡単に言えば、あなたはそれを単に次のように考えることができます データサイエンス 生物学のために。

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

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

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

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

ヒトゲノムの最新のGFF3ファイルは以下からダウンロードできます。 ここに 。ザ・ README このディレクトリにあるファイルは、このデータ形式の簡単な説明を提供し、より完全な仕様が見つかります ここに 。

我々は使用するだろう パンダ 、GFF3ファイルを操作および理解するための、高速で柔軟で表現力豊かなデータ構造を提供するSciPyスタックの主要コンポーネント。

セットアップ

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

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

-b bash行のフラグは、バッチモードで実行するように指示します。上記のコマンドを使用して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つ

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 、 ソース 、 タイプ 、 開始 、 終わり 、 スコア 、 ビーチ 、 段階 、 属性 。それらのいくつかは非常に理解しやすいです。例として最初の行を取り上げます。

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

これは、1番目の塩基から24,895,622番目の塩基まで始まる1のseqidを持つ最初の染色体の注釈です。

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

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

最後の属性列には、1番染色体にもCM000663.2、chr1、NC_000001.11の3つのエイリアス名があることが示されています。これは基本的にGFF3ファイルの外観ですが、1行ずつ検査することはないので、ファイル全体を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=' ', comment='#', low_memory=False, header=None, names=col_names)

上記の最後の行は、GFF3ファイル全体にpandas.read_csvをロードします。方法。

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

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

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

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

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

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

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

パンダにもSeriesがありますおよびPanelそれぞれ1Dおよび3Dデータのデータ構造。を参照してください ドキュメンテーション パンダのデータ構造の紹介。

.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)を使用して、長い文字列を省略しないようにパンダを設定できます。さらに、パンダにはたくさんあります オプション カスタマイズすることができます。

次に、.infoを使用してこのデータフレームに関する基本情報を取得しましょう。方法。

In [20]: df.info() 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つの方法です。別の方法は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'から評価された値を調べると、それは一連のTrueです。およびFalse dfと同じインデックスを持つ各エントリの値。 df[]に渡します対応する値がTrueであるエントリのみを返します。

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

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

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

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

次に、各エントリの長さを計算し、それを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.+?);') 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.+?); 、+? +の代わりに使用されます欲張りではなく、最初のセミコロンで検索を停止させたいからです。それ以外の場合、結果は最後のセミコロンまで一致します。

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

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

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

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

RE_GENE_ID = re.compile(r'gene_id=(?PENSG.+?);') 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.+?);') 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遺伝子を詳しく見ると、実際にはすべて異なっていることがわかります。

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

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

次のどれが卸売りだけに影響を与える5つの競争要因の1つではありません

ここでのポイントは、遺伝子名がすべての遺伝子IDに固有であるとは限らず、複数の遺伝子によって共有される名前の約0.15％であることを知っていることです。

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

ゲノムの不完全性を理解しようとしたときと同じように、lengthを簡単に追加できます。列から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()長さの値に基づいていくつかの簡単な統計を計算します。

• 遺伝子の平均の長さは約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ファイル内にありますが、この投稿の範囲外です。

染色体間の遺伝子分布

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

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）には遺伝子がないように見えることに注意してください。これは真実ではありません。

最初のDataFrameでもう少しフィルタリング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...

この例は、フィルタリング中に2つの条件を&で組み合わせる方法も示しています。 「または」の論理演算子は|になります。

各条件を括弧で囲む必要があり、Pandasの構文のこの部分は、リテラルであるPythonとは異なることに注意してくださいandおよびor。

次に、gdfを借りましょう各染色体の長さのソースとしての前のセクションのDataFrame：

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

明確にするために、分析に関係のない列は削除されています。

はい、.drop列のリストを取得して、1回の操作でそれらをまとめて削除することもできます。

MTのseqidを持つ行がまだそこにあることに注意してください。後で戻ります。次に実行する操作は、seqidの値に基づいて2つのデータセットをマージすることです。

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

chr_gene_counts以降はまだSeriesオブジェクトであり、マージ操作をサポートしていません。最初に.to_frameを使用してDataFrameオブジェクトに変換する必要があります。

.reset_index()元のインデックス（つまり、seqid）を新しい列に変換し、現在のインデックスを0ベースの増分番号としてリセットします。

cdf.head(2)からの出力それがどのように見えるかを示しています。次に、.mergeを使用しましたseqid列（on='seqid'）で2つのDataFrameを結合するメソッド。

マージ後gdfおよびcdf、MTエントリがまだありません。これは、デフォルトで.mergeが原因です。内部結合を操作しますが、左結合、右結合、または外部結合は、howを調整することで使用できます。パラメータ。

を参照してください ドキュメンテーション 詳細については。

後で、関連する.joinもあることに気付くかもしれません。方法。 .mergeおよび.join似ていますが、APIが異なります。

関係者によると ドキュメンテーション 言う

関連するDataFrame.joinメソッドは、index-on-indexおよびindex-on-column（s）結合に内部的にマージを使用しますが、共通の列で結合しようとするのではなく、デフォルトでインデックスで結合します（マージのデフォルトの動作）。インデックスで結合している場合は、DataFrame.joinを使用して入力を節約することをお勧めします。

基本的に、.mergeより汎用的で、.joinによって使用されます。

最後に、染色体length間の相関を計算する準備ができました。および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

デフォルトでは.corrを計算します ピアソン相関 データフレーム内の列のすべてのペア間。

ただし、この場合、列のペアは1つしかなく、相関は正の0.73であることがわかります。

言い換えれば、染色体が大きいほど、より多くの遺伝子を持っている可能性が高くなります。値のペアをlengthで並べ替えた後、2つの列もプロットしてみましょう。

ax = merged[['length', 'gene_count']].sort_values('length').plot(x='length', y='gene_count', style='o-') # 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))

上の画像に見られるように、それは全体的に正の相関関係ですが、すべての染色体に当てはまるわけではありません。特に、17番染色体、16番染色体、15番染色体、14番染色体、13番染色体の場合、相関は実際には負であり、染色体のサイズが大きくなるにつれて染色体上の遺伝子の数が減少することを意味します。

調査結果と将来の研究

これで、SciPyスタックを使用したGFF3形式のヒトゲノムのアノテーションファイルの操作に関するチュートリアルは終了です。主に使用しているツールには、IPython、Pandas、matplotlibなどがあります。チュートリアルでは、パンダで最も一般的で便利な操作のいくつかを学んだだけでなく、ゲノムに関するいくつかの非常に興味深い質問にも答えました。要約すれば：

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

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

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