Juliaでバイオインフォマティクス（ゲノムデータ編）

この記事は Julia Advent Calendar 2018 の7日目の記事です。

TL; DL

JuliaでバイオインフォマティクスができるBioJuliaについて紹介します。

ゲノムデータって何？

ゲノムって何？

雑にまとめると、DNAという物質でできた、生命の設計図です。

データとしては、A/C/G/Tの４種類のアルファベットからなる文字列で、この文字列が染色体ごとにあります。
「染色体を表す文字列（染色体名; chr1など）」と「その染色体に対応する文字列の何文字目かを表す整数」によってゲノム上の位置や区間を定義することができます。

どんなデータがあるの？

• 配列データ
• ゲノム座標上の区間・シグナル
• アラインメント/マッピング（２本の文字列間でどの文字とどの文字が対応しているかを表す）
• 変異（参照ゲノム配列を基準にした際の、個人の持つゲノム配列における差分）

BioJuliaについて

BioJuliaは生物学的データを取り扱うパッケージを開発するJuliaのコニュニティ。Juliaでお馴染みの@bicycle1885さんも古くからコミットしています。

バイオ系は様々な種類のデータを表現するために様々な種類のフォーマットが使われます。BioJuliaはそれらの生物系データの主なフォーマットのI/Oを提供しています。
また、ゲノムデータでよく使われる計算（例: 相補鎖を作る、配列アラインメント）も実装されています。

どんなパッケージがあるの？

ゲノムデータを扱うパッケージとして主なものを挙げます。フォーマットの仕様は検索すれば出てくるので、説明を省きます。

データ	主なファイルフォーマット	BioJulia のパッケージ
生物配列	FASTA, FASTQ, 2bit	BioSequences
ゲノム上の特徴	BED, BigWig, BigBed, GFF3	GenomicFeatures
配列アラインメント	SAM, BAM	BioAlignments
遺伝的変異	VCF, BCF	GeneticVariation

BioJuliaのページ には他にも様々なパッケージがあります。

例

BioSequences

julia> using BioSequences

# 文字列の前に `dna` を前につければDNA配列が作れる
julia> seq = dna"CACTGT"
6nt DNA Sequence:
CACTGTT

julia> typeof(seq)
BioSequence{DNAAlphabet{4}}

# 相補鎖（AとT、GとCが相補的というやつ）を計算できる
julia> reverse_complement(seq)
6nt DNA Sequence:
ACAGTG

セントラルドグマ（DNA-(転写)->RNA-(翻訳)->タンパク質）の再現もできます。

# DNA→RNA
julia> rna = convert(RNASequence, seq)
6nt RNA Sequence:
CACUGU

# RNA→タンパク質
julia> BioSequences.translate(rna)
2aa Amino Acid Sequence:
HC

# DNAからも直接翻訳したアミノ酸配列が得られます
julia> BioSequences.translate(seq)
2aa Amino Acid Sequence:
HC

その他に、マッチする配列を検索するといったこともできます。

GenomicFeatures

東北メディカルメガバンクが所収する、血球細胞のクロマチンアクセシビリティのデータを可視化してみます。クロマチンアクセシビリティは、ゲノムの特定の領域において分子装置などのアクセスしやすさで、臓器・細胞における遺伝子のON/OFFを反映します（クロマチンアクセシビリティが高い遺伝子はONであるｋじょとが多い）。

BigWigフォーマットで格納されています。これは、ゲノム座標上のシグナルをRun length encodingで記述したものです。データはこちらからダウンロードできます。

試しに、T細胞の表面抗原であるCD4遺伝子の周辺領域でのシグナルをパースし、プロットしてみました。

julia> using GenomicFeatures

julia> function readBigWigValues(pathbw::String, interval::Interval)
           open(BigWig.Reader, pathbw) do reader
               return(BigWig.values(reader, interval))
           end
       end
readBigWigValues (generic function with 1 method)

julia> pathbw = "tommo-genome_accessibility-20181105-259pe_mapq00-MERGED.bigWig"
"Downloads/tommo-genome_accessibility-20181105-259pe_mapq00-MERGED.bigWig";

julia> chrom = "12";

julia> leftpos = 6783678;

julia> rightpos = 6821003;

julia> y = readBigWigValues(pathbw, Interval(chrom, leftpos, rightpos))
37326-element Array{Float32,1}:
 23.0
 23.0
 23.0
 23.0
 23.0
 23.0
 23.0
 23.0
 23.0
 23.0
  ⋮
 12.0
 12.0
 12.0
 12.0
 12.0
 12.0
 12.0
 12.0
 12.0
 12.0

julia> using Plots

julia> plot(leftpos:rightpos,y)

図の右側のピークが立っている部分は、CD4遺伝子の転写開始点に相当し、クロマチンアクセシビリティが高くなっていることがわかります。

その他に、２つの区間集合の集合演算などもできます。

BioAlignments

今度はBAMファイルの可視化をしてみます。BAMファイルはゲノムのどの位置と対応するか（マッピング、アラインメント）を表します。ここでは、ゲノム座標ごとにいくつの配列がアラインメントされているかのカウントを表示します。

データは、１細胞RNA-seqの一種RamDA-seqによるマウス胚性幹細胞のデータから、19番染色体（短くて扱いやすい）だけにして扱いやすくしたテスト用データです。

ここではRNAスプライシングによる不連続なアラインメント（スプリットアラインメント）を考慮し、スプリットが起こっているところでは配列をカウントしないようにします。そのために自分で関数を書いています。

julia> using BioAlignments

julia> function bamToCoverage_cigarAware_base(bamReader::BAM.Reader, chrom::String, leftpos::Int64, rightpos::Int64)
    cov = zeros(Int, rightpos - leftpos + 1)
    for record in eachoverlap(bamReader, chrom, leftpos:rightpos)
        # skip unmapped records
        if ! BAM.ismapped(record)
            continue
        end

        readLeftPos = BAM.position(record)
        cigarRle = BAM.cigar_rle(record)
        offset = 0
        for i in 1:length(cigarRle[1])
            if cigarRle[1][i] == OP_MATCH && leftpos ≤ readLeftPos + offset ≤ rightpos
                for j in 1:min(cigarRle[2][i], rightpos - readLeftPos - offset + 1)
                    cov[readLeftPos + offset + j - leftpos] += 1
                end
            end
            offset += cigarRle[2][i]
        end
    end
    return(cov)
end

julia> function bamToCoverage_cigarAware(pathBam::String, chrom::String, leftpos::Int64, rightpos::Int64)
    open(BAM.Reader, pathBam, index=pathBam*".bai") do reader
        return(bamToCoverage_cigarAware_base(reader, chrom, leftpos, rightpos))
    end
end


julia> pathBamFile = "RamDA_00h_A04.uniq.q40.chr19.bam";

julia> chrom = "chr19";

julia> leftpos = 3205801;

julia> rightpos = 3206000;

julia> cov = bamToCoverage_cigarAware(pathBamFile, chrom, leftpos, rightpos)
200-element Array{Int64,1}:
 1
 1
 1
 1
 1
 1
 1
 1
 1
 1
 ⋮
 0
 0
 0
 0
 0
 0
 0
 0
 0
 0

julia> Plots.plot(leftpos:rightpos, cov)

BioAlignmentsはこの他に、アラインメント自体を計算する関数も備えています。

おわりに

いまのところ足りないもの

最後にいまのところ足りてないと思うものを挙げます。ここに足りないものは、まずゲノムブラウザですよね。
ゲノムブラウザはゲノム座標に沿って、様々な種類のゲノムデータを並行させて表示するビューワーです。

Dalliance というゲノムブラウザをJuliaから呼び出せるGenomeBrowsers.jlというパッケージがあるのですが、現状インストールできず、またソースを見た感じインストールできたとしても動きません（Julia 1.0に達する前の燃え尽き症候群な気がします）。

また、2016年にはJBrowseなどのゲノムブラウザとBioJuliaのパッケージ群との連携の可能性が述べられていましたが、現状どうなっているか不明です。

結論

元々@bicycle1885さんにBioJuliaに薦めてもらって、最終的にBioJuliaを活用して論文を書くことができました（論文、Julia Computingの紹介記事）。使い勝手がよく今後も使っていきたいパッケージ群です。

ゲノムデータをいじる予定がある方は、自分でパーサー書く前に、ぜひBioJuliaを試してみてほしいです。