Python入門

Pythonとは

Pythonは、1991年にGuido van Rossum(グイド・ヴァンロッサム)氏により開発されたプログラミング言語です。C言語などのように機械に近い言語を「低級言語」と呼びますが、Pythonは比較的人間が理解しやすい「高級言語」に分類されます。Pythonは高級言語の中でも、特に習得が容易であるため、プログラミング初心者に向いている言語であると言われています。しかしながら、初心者だけが使用するプログラミング言語というわけではありません。Pythonは多くの企業で採用されており、実践的なプログラミング言語でもあります。特にデータ分析やWEBアプリ開発の分野で、よく利用されています。

Pythonを選ぶ理由

学習コストが低い

Pythonは比較的シンプルな文法ですので、学習にかかる時間も少なく済みます。

汎用的である

Pythonは他の言語と比べると多くの用途に利用できます。データ分析の分野に興味があるのであれば、Pythonが第一選択肢になるかと思います。WEBアプリ開発の分野でも多く利用されています。よく比較される言語であるPerlやRubyやPHPの場合、WEBアプリ開発には向いていますが、データ分析を本格的に行うには向いていません。

WEB上に情報が多い

これは、Pythonのユーザー数が多いことに起因します。ユーザー数が多ければ、その分、情報を発信する人が多くなります。

ライブラリが充実

ライブラリとは、よく使用される機能をひとまとまりにしたプログラムです。ライブラリが充実していれば、開発時間を短縮することが可能です。これも、Pythonのユーザー数が多いことに起因します。

バイオインフォマティクス解析でよく使う(かもしれない)gitコマンドのチートシート

現在いるブランチ名を表示する

$ git branch --contains

変更されたファイルすべてをaddする

$ git add -A

変更内容をコミット

$ git commit -m "comment"

masterブランチに移動

$ git checkout master

branch1ブランチを現在のブランチにマージ

$ git merge branch1

ローカルのmasterブランチをリモートのmasterにpush

$ git push origin master

ローカルのbranch1ブランチをリモートのmasterにpush

$ git push origin branch1:master

登録されているリモートリポジトリの確認

$ git remote -v

バイオインフォマティクス解析でよく使う(かもしれない)awkコマンドのチートシート

2列目の値を表示する

$ awk '{print $2}’ [ファイル名]

2列目の値が”foo”に一致する行だけを表示する

$ awk '$2 == "foo"' [ファイル名]

2列目の値の合計を表示する

$ awk '{a+=$2} END{print a;}' [ファイル名]

STARとhtseq-countでRNAseq解析を行う

データの準備

解析対象とするfastqファイルをダウンロードしていきます。
SRA Toolkitをインストールしていない場合には、以下のコマンドでインストールできます。

$ wget https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/current/sratoolkit.current-ubuntu64.tar.gz
$ tar zxvf sratoolkit.current-ubuntu64.tar.gz

binフォルダにパスを通しておいてください。
以下のコマンドでSRR3485766のfastqファイルをダウンロードできます。

$ fastq-dump -A SRR3485766 --split-files 

これで、カレントディレクトリに SRR3485766_1.fastqとSRR3485766_2.fastqが作成されました。

マッピングの準備

参照配列をダウンロードします。

$ wget ftp://ftp.ensembl.org/pub/release-94/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
$ gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
$ wget ftp://ftp.ensembl.org/pub/release-94/gtf/homo_sapiens/Homo_sapiens.GRCh38.94.gtf.gz
$ gunzip Homo_sapiens.GRCh38.94.gtf.gz

次に、STARのインデックスを作成していきます。
STARをインストールしていない場合には、以下のコマンドでインストールできます。

$ wget https://github.com/alexdobin/STAR/archive/2.6.1d.tar.gz
$ tar zxvf 2.6.1d.tar.gz

bin/Linux_x86_64にパスを通しておきます。
以下のコマンドでインデックスを作成します。

$ STAR \
  --runThreadN 8 \
  --runMode genomeGenerate \
  --genomeDir genome \
  --genomeFastaFiles Homo_sapiens.GRCh38.dna.primary_assembly.fa \
  --sjdbGTFfile Homo_sapiens.GRCh38.94.gtf 

マッピング

以下のコマンドで、マッピングをします。

$ STAR --genomeDir genome/ \
  --runThreadN 8 \
  --readFilesIn SRR3485766_1.fastq SRR3485766_2.fastq \
  --outSAMtype BAM SortedByCoordinate \
  --outFileNamePrefix sample1

リードカウント

HTSeqをインストールしていない場合には、以下のコマンドでインストールできます。

$ pip install numpy
$ pip install htseq

以下のコマンドで、リードカウントをします。

$ htseq-count -f bam -r pos -t exon \
  sample1Aligned.sortedByCoord.out.bam \
  Homo_sapiens.GRCh38.94.gtf \
  > result.txt

result.txtにリードカウント結果が得られました。

中国で「ゲノム編集の双子」が誕生した件についてのまとめ

2018年11月26日、中国の研究者 賀建奎副教授がゲノム編集の技術を人の受精卵に応用し、双子の女児を誕生させたと主張して話題になっています。

問題の動画:

HIVに耐性を持った女児が生まれたと主張

遺伝子編集ベイビー

ヒトの受精卵をゲノム編集で操作し子供を誕生させることは、倫理、技術の両面から問題が指摘されており、今までに誕生した例はありませんでした。今回の主張が事実だとすると、世界初の誕生例となります。

賀建奎副教授について

賀建奎副教授が所属する南方科技大学のWEBページの紹介によると、

  • 2006年、中国科学技術大学 現代物理学 学士号取得
  • 2010年、ライス大学 生物物理学 博士号取得
  • 2011-2012年 スタンフォード大学 ポスドク

という経歴です。その後、南方科技大学に就職したとされています。
また、賀建奎副教授はビジネスにも積極的で、7社の株主、6社の代表取締役を務めているようです。その中でも、最初に設立したDirect Genomics社は2018年4月に2.18億元を調達するなど、注目を集めていました(科学网)。

批判

南方科技大学は声明を発表し、賀建奎副教授は現在休職中(2018年2月—2021年1月)であり、今回の件について大学は関与していないことを説明しました。今後、独立委員会を設置し、調査を進める予定であるとのことです。
また、中国の科学者の間からも批判の声が上がっており、122人が共同声明を発表し、「人体実験」を行ったことに対して反対を表明しました。

Second International Summit on Human Genome Editingでの発表

2018年11月28日、香港大学で開かれた国際ゲノム編集会議で賀建奎副教授が発表を行いました。この発表は、今回の騒動以前から予定されていたものであり、元の発表内容を変更して行われました。発表では、マウスにおけるCCR5ノックアウトについての説明、そして今回の件のヒトにおけるCCR5ノックアウトについての説明がありました。その後、参加者から質問を受け付けました。その中で、遺伝病を持った患者たちを助けるために、今回の成果が役に立つことを強調しました。また、今回生まれた双子以外にも、もう一人妊娠している女性がいることを明らかにしました。

参考:贺建奎的“基因编辑婴儿”是一项突破吗?不,可能是倒退

バイオインフォマティクスにおけるパイプラインの構築にはどのフレームワークを使用すれば良いか

バイオインフォマティクス解析においては、多数のOSS(オープンソースソフトウェア)を組み合わせてパイプラインを構築することがよくあります。パイプラインの構築は、自分の好きな言語を選んで自己流で行うことも可能ですが、パイプライン構築用のフレームワークを使用することで、可読性が高く、簡潔なプログラムを書くことができるようになります。
多くのフレームワークではDockerコンテナに対応しているため、多数のOSSの実行環境を整備するといった手間のかかる作業を短縮することができます。また、Dockerコンテナを使用すると、パイプラインのポータビリティが高まりますので、他の研究者へパイプラインを渡したり、クラウド上でパイプラインを実行したりすることも容易になります(再現性の向上)。
パイプライン構築用のフレームワークは、awesome-pipelineにまとめられています。たくさんありますので、そのうち主なものを以下にまとめました。

snakemake

Pythonベースのフレームワーク。Pythonの記法が使えるので、柔軟なワークフローが設計できる。ただ、柔軟なワークフローが設計できる分、タスクごとにDockerで実行環境を分離する考え方と相性が悪く、パイプラインのポータビリティはあまり良くない印象。

例)

rule targets:
    input:
        "plots/dataset1.pdf",
        "plots/dataset2.pdf"

rule plot:
    input:
        "raw/{dataset}.csv"
    output:
        "plots/{dataset}.pdf"
    shell:
        "somecommand {input} {output}"

Common Workflow Language(CWL)

ワークフローを記述するための統一言語。CWLを実行するためのソフトウェアはcwltoolArvadostoil等たくさん開発されている。

例)

#!/usr/bin/env cwl-runner

cwlVersion: v1.0
class: CommandLineTool
baseCommand: echo
inputs:
  message:
    type: string
    inputBinding:
      position: 1
outputs: []

Nextflow

CWLと目指す方向はほぼ一緒であるが、CWLは言語仕様と実行エンジンが分離しているのに対して、Nextflowは同一である。その分、Nextflowの方がCWLと比べて柔軟にワークフローを記述できる印象。

例)

params.query = "$HOME/sample.fa"
params.db = "$HOME/tools/blast-db/pdb/pdb"

process blast {
    output:
     file top_hits

    """
    blastp -query ${params.query} -db ${params.db} -outfmt 6 \
    | head -n 10 \
    | cut -f 2 > top_hits
    """
}

process extract {
    input:
     file top_hits
    output:
     file sequences

    "blastdbcmd -db ${params.db} -entry_batch $top_hits > sequences"
}

process align {
    input:
     file sequences
    echo true

    "t_coffee $sequences 2>&- | tee align_result"
}

メモ

NextflowのGitHubスター数が635、CWLのGitHubスター数が815(2018年11月19日)であることから、CWLがやや優勢?

バイオインフォマティクスを始めるのにどの言語を勉強したら良いか

バイオインフォマティクスをこれから始める人にどの言語を使用したら良いかと質問されることがよくある。周りのバイオインフォマティシャンに質問してみると、たぶん返ってくる答えはPython, Perl, Rubyのいずれかだと思う。結論からいうと、2018年現在でバイオインフォマティクスをこれから始める人にふさわしい言語はPython一択であると考えられる。

Pythonはユーザー数が多い

それぞれの言語のユーザー数がどれくらいいるか正確に把握することは難しいが、参考までにGoogleトレンドの結果を示す。


全世界:

日本:

世界的には2011年後頃からPerlを逆転してPythonがトップになっている。日本でも2014年頃にはPythonがPerlを抜いている。

なぜユーザー数が多い言語を選ぶべきか

Pythonのユーザー数が多いことは理解いただけたかと思うが、なぜユーザー数が多い言語を選ぶべきかを説明する。

理由1. ドキュメントが豊富

ドキュメントが豊富である一例として、プログラミングに関するQ&Aサイトである「Stack Overflow」の検索結果を示す。

Python 1,059,777 件
Perl 60,952件
Ruby  199,959 件

圧倒的にPythonの件数が多くなっている。WEB上に情報が多ければ、問題が起きた時にググるだけで解決する可能性が高い。

理由2. ライブラリーが充実

ライブラリーが充実してる一例として、バイオインフォマティクス向けのライブラリーであるBioPython, BioPerl, BioRubyを比較してみる。まずはそのライブラリーの開発に貢献した人数である「Contributer」の数を見てみる。

BioPython 208人
BioPerl 62人
BioRuby 30人

BioPythonが最も多くの人が関わっている。次に、直近1年の更新頻度を見てみる。

BioPython:

BioPerl:

BioRuby:

縦軸に注意してもらいたい。圧倒的にBioPythonの更新頻度が多くなっている。
ユーザー数の多い言語のライブラリーほど、多くの人が関わり、更新頻度が高くなるので、機能が充実している&バグが少ない可能性が高い。

まとめ

以上のように、バイオインフォマティクスをこれから始める人にふさわしい言語はユーザー数の多いPythonであると考えられる。過去に開発されたプログラムはPerlだからPerlを使うべきという人がいるが、進歩の早いこの分野においては、10年前のプログラムは使い物にならないことが多い。また、Rubyを勧める人の言い分で「Rubyは日本語ドキュメントが充実している」という人がいるが、主にWEB開発に関するドキュメントであり、バイオインフォマティクスにはほとんど役に立たない。今後他の言語が主流になる可能性は否定しないが、現状としてはPythonをお勧めする。

ロングリードでのアセンブルは遺伝子予測結果にエラーが多く含まれる

ゲノム配列がまだ決定されていない生物種について新規にゲノム配列を決定する「de novoアセンブリ解析」では、第2世代シーケンサー(Illumin等)と比較して長いリード長が得られる第3世代シーケンサー(PacBio, NanoPore等)を使用することが多い。しかしながら、第3世代シーケンサーではリードにエラーが多く含まれており、そのリードをもとにアセンブルすると、遺伝子領域に多くのギャップが含まれた結果が得られてしまう。第2世代シーケンサーの結果を使用してアセンブル結果のエラー補正を実施した場合でも、多くのギャップが残っている。

リファレンス Gapが含まれる遺伝子の数 Gapが含まれる遺伝子の割合
GRCh38(コントロール) 161
Nanopore (Jain et al) 7022 33 %
PacBio1 (Pendleton et al) 16630 80 %
PacBio2 (Koren et al)  5272 25 %

参考情報

Mind the gaps – ignoring errors in long read assemblies critically affects protein prediction

メタゲノム解析

メタゲノム解析とは

メタゲノム解析とは、微生物群を直接シーケンスし、網羅的に菌種や遺伝子を特定する解析のことです。糞便や皮膚、海洋や土壌等といったサンプルを対象として研究が進められています。解析の手法には大きく分けて、全ゲノムを対象としてシーケンスする「ショットガンメタゲノム解析」と16S rRNAのみを対象としてシーケンスする「16S rRNAメタゲノム解析」の2種類があります。現在は主にショートリードシーケンサーが使用されていますが、今後はロングリードシーケンサーを使用することで、より精度の高い解析が可能になることが期待されます。

ショットガンメタゲノム解析

ショットガンメタゲノム解析におけるバイオインフォマティクス解析の手法は、大きく分けて2つあります。1つはリファレンス配列を使用せずアセンブリを行う手法で、「MetaSPAdes」や「MEGAHIT」といったソフトウェアがあります。もう1つは、リファレンス配列に対してマッピングをする手法で、「MetaPhlAn」といったソフトウェアがあります。ヒトの糞便サンプル等といった研究が比較的進んでいる分野については、マッピングの手法が使われることが多いようです。

16S rRNAメタゲノム解析

16S rRNAメタゲノム解析は16S rRNAのみをシーケンス対象とすることで、解像度は落ちる分、低コストで解析を実施することが可能です。16S rRNAメタゲノム解析には「QIIME」といったソフトウェアがよく使われます。

バイオインフォマティクス解析でよく使われるファイル形式

バイオインフォマティクス解析では様々なファイル形式が使用されています。ここでは、主要なファイル形式を紹介します。