需求

前幾天談到的東西幾乎都是與序列比對有關，特別是BLAST之類的比對方法及結果。今天要談的是另外一些比對工具產生的結果，以及怎麼使用Julia底下的工具來處理它們。

不同的檔案格式源自於不同的需求

在講今天想要分享的工具前，我必須先大概地介紹一下這個工具能處理的那些不同資料的格式，詳細資料就讓有興趣的朋友們自行去搜尋了：

• Fastq: 這種格式的檔案我在前幾天[從Fastq出發到abundance table]一文裡頭有介紹過，存放的是序列及序列對應的quality score。
• SAM/BAM：SAM這種格式的文字檔存放了將Fastq格式的檔案比對到特定的參考序列後的相關資訊，像是哪條序列比對到參考序列的哪裡、比對上的長度、比對的分數、比對狀態標籤以及其他細節；BAM格式的檔案則是將SAM格式的檔案轉換為二進位的格式儲存。
• BED：這種格式的檔案主要存放了參考序列中每段區間所對應的染色體、區間座標、正/反股標記以及這段區間對應的實質生物意義，可以用來做註解也可用來取參考序列中區間內的序列片段。
• GTF：這個格式的檔案與BED格式的檔案功能很接近，但主要是用來做註解使用，其中擺放的註解資訊通常也比BED檔都很多。
• VCF:這種格式的檔案主要是用來標示最原始的Fastq檔案經比對、註解以及變異位點辨識等過程後的結果。

OpenGene Project

OpenGene.jl這個Julia的package正好提供了我們相當方便的工具來處理上述提到的各種檔案，我在下面分別舉些例子：

• 安裝及基本使用

Pkg.add("OpenGene")
using OpenGene
# 產生一段DNA序列
seq = dna("AAATTTCCCGGGATCGATCGATCG")
# 取得這段序列的倒序列
~seq
# 轉錄這段序列
transcribe(seq)

• 讀寫fastq/fasta檔

using OpenGene
istream = fastq_open("input.fastq.gz")
ostream = fastq_open("output.fastq.gz","w")

while (fq = fastq_read(istream))!=false
    fastq_write(ostream, fq)
end

close(ostream)

• 讀/寫Illumina定序儀產生的成對Fastq檔

using OpenGene

istream = fastq_open_pair("R1.fastq.gz", "R2.fastq.gz")
ostream = fastq_open_pair("Out.R1.fastq.gz","Out.R2.fastq.gz","w")

while (pair = fastq_read_pair(istream))!=false
    fastq_write_pair(ostream, pair)
end

close(ostream)

• 讀/寫BED檔

using OpenGene
intervals = bed_read_intervals("in.bed")
bed_write_intervals("out.bed",intervals)

• 讀VCF檔，並對讀進來的資料進行處理再寫入新的VCF檔

using OpenGene

v1 = vcf_read("in_1.vcf")
v2 = vcf_read("in_2.vcf")
v1 = vcf_read("v1.vcf")
v2 = vcf_read("v2.vcf")

# 按位置merge
v_merge = v1 + v2

# 按位置取交集
v_intersect = v1 * v2

# 移除掉v1中同時也出現在v2裡的紀錄
v_minus = v1 - v2

# 寫入vcf檔
vcf_write("minus.vcf", v_minus)

• 讀/寫GTF檔

using OpenGene

gtfobj = gtf_read("in.gtf")
gtf_write("out.gtf", gtfobj)
gtfobj, stream = gtf_read("in.gtf", loaddata = false)
while (row = gtf_read_row(stream)) != false
    # do something with row ...
end

以上大概就是這個Package的功能，不過僅僅只是相當粗略的介紹，還待我花更多的時間去摸索再與各位分享了。