Python: OpenCVを使用して顔検出をする

OpenCVを使うと、簡単に顔を検出することができます。

import cv2

cascade = cv2.CascadeClassifier('/usr/local/lib64/python3.7/site-packages/cv2/data/haarcascade_frontalface_default.xml')
img = cv2.imread('img/img.jpg')
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
rects = cascade.detectMultiScale(gray, 1.1, 3,0, (20,20))

if len(rects) > 0:
    for rect in rects:
        cv2.rectangle(img, tuple(rect[0:2]), tuple(rect[0:2]+rect[2:4]), (0, 0,255), thickness=2)
else:
    print("検出なし")

cv2.imwrite('img/img_cascade.jpg', img)

分類器は、

cascade = cv2.CascadeClassifier('/usr/local/lib64/python3.7/site-packages/cv2/data/haarcascade_frontalface_default.xml')

として指定していますが、分類器のファイルの場所は以下のようにして確認できます。

import cv2
print(cv2.__path__)
#['/usr/local/lib64/python3.7/site-packages/cv2']

フォルダの中を確認してみますと、たくさんの分類器があることがわかります。今回はそのうち、「haarcascade_frontalface_default.xml」を使用していきます。

cd /usr/local/lib64/python3.7/site-packages/cv2/data
ls
#__init__.py					haarcascade_frontalface_alt2.xml		haarcascade_profileface.xml
#__pycache__					haarcascade_frontalface_alt_tree.xml		haarcascade_righteye_2splits.xml
#haarcascade_eye.xml				haarcascade_frontalface_default.xml		haarcascade_russian_plate_number.xml
#haarcascade_eye_tree_eyeglasses.xml		haarcascade_fullbody.xml			haarcascade_smile.xml
#haarcascade_frontalcatface.xml			haarcascade_lefteye_2splits.xml			haarcascade_upperbody.xml
#haarcascade_frontalcatface_extended.xml		haarcascade_licence_plate_rus_16stages.xml
#haarcascade_frontalface_alt.xml			haarcascade_lowerbody.xml

分類器のパスを間違えていると、以下のようなエラーが出ます。

Traceback (most recent call last):
  File "[*.py]", line 6, in 
    rects = cascade.detectMultiScale(gray, 1.1, 3,0, (20,20))
cv2.error: OpenCV(4.0.0) /io/opencv/modules/objdetect/src/cascadedetect.cpp:1658: error: (-215:Assertion failed) !empty() in function 'detectMultiScale'

顔検出結果

顔検出の結果は以下の通りです。画像は写真素材 足成からいただきました。

元画像1

顔検出結果1

元画像2

顔検出結果2

元画像3

顔検出結果3

モザイク処理をする

顔検出した箇所についてぼかしを入れることもできます。

import cv2

cascade = cv2.CascadeClassifier('/usr/local/lib64/python3.7/site-packages/cv2/data/haarcascade_frontalface_default.xml')
img = cv2.imread('img/img.jpg')
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
rects = cascade.detectMultiScale(gray, 1.1, 3,0, (20,20))

if len(rects) > 0:
    for rect in rects:
        x, y, w, h = rect
	face = img[y:y+h, x:x+w]
        dst = cv2.GaussianBlur(face, (25, 25), 0)
        img[y:y+h, x:x+w] = dst
else:
    print("検出なし")

cv2.imwrite('img/img_mosaic.jpg', img)

モザイク処理結果1

モザイク処理結果2

モザイク処理結果3

Python: OpenCVを使用してエッジを検出する

画像からエッジを検出する代表的な手法として、Cannyエッジ検出というものがあります。OpenCVを使用すると、簡単にCannyエッジ検出を行うことができます。

cv2.Canny(image, threshold1, threshold2)

Cannyエッジ検出では、画像データに加えて、2つの引数threshold1とthreshold2が必須です。この2つの値はエッジ検出の際に以下のように使用されます。

threshold2より大きい場合: エッジとして検出される
threshold1より大きくthreshold2より小さい場合: エッジと隣接している場合には、エッジとして検出される
threshold1より小さい場合: エッジとして検出されない

エッジ検出結果

元画像

threshold1=100, threshold2=100

import cv2

img = cv2.imread('img/img.jpg')
dst = cv2.Canny(img, 100, 100)
cv2.imwrite('img/img_canny.jpg', dst)

結果

少し線が多すぎる印象です。

threshold1=100, threshold2=300

import cv2

img = cv2.imread('img/img.jpg')
dst = cv2.Canny(img, 100, 300)
cv2.imwrite('img/img_canny.jpg', dst)

結果

いい感じにエッジが検出できました。

threshold1=100, threshold2=600

import cv2

img = cv2.imread('img/img.jpg')
dst = cv2.Canny(img, 100, 600)
cv2.imwrite('img/img_canny.jpg', dst)

結果

少し線が少ないでしょうか。

threshold1=600, threshold2=600

import cv2

img = cv2.imread('img/img.jpg')
dst = cv2.Canny(img, 600, 600)
cv2.imwrite('img/img_canny.jpg', dst)

結果

threshold1を上げると、線のつながりがなくなってしまいました。

Python: OpenCVのGaussianBlurを使用して画像にぼかしをかける

ガウシアンフィルタは、画像の平滑化に使われるフィルタの1つです。ガウス分布を利用して、注目画素からの距離に応じて近傍の画素値に重みをかけます。

インストール

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

pip install opencv-python

ぼかしをかける

元画像

カーネルサイズの調整

カーネル=(1, 1)、標準偏差3

import cv2

img = cv2.imread('img/img.jpg')
dst = cv2.GaussianBlur(img, (1, 1), 3)
cv2.imwrite('img/img_gauss.jpg', dst)

結果

カーネル=(5, 5)、標準偏差3

import cv2

img = cv2.imread('img/img.jpg')
dst = cv2.GaussianBlur(img, (5, 5), 3)
cv2.imwrite('img/img_gauss.jpg', dst)

結果

カーネル=(15, 15)、標準偏差3

import cv2

img = cv2.imread('img/img.jpg')
dst = cv2.GaussianBlur(img, (15, 15), 3)
cv2.imwrite('img/img_gauss.jpg', dst)

結果

カーネルサイズを大きくすると、ぼかしの具合が強くなります。

標準偏差の調整

カーネル=(15, 15)、標準偏差1

import cv2

img = cv2.imread('img/img.jpg')
dst = cv2.GaussianBlur(img, (15, 15), 1)
cv2.imwrite('img/img_gauss.jpg', dst)

結果

カーネル=(15, 15)、標準偏差2

import cv2

img = cv2.imread('img/img.jpg')
dst = cv2.GaussianBlur(img, (15, 15), 2)
cv2.imwrite('img/img_gauss.jpg', dst)

結果

カーネル=(15, 15)、標準偏差3

import cv2

img = cv2.imread('img/img.jpg')
dst = cv2.GaussianBlur(img, (15, 15), 3)
cv2.imwrite('img/img_gauss.jpg', dst)

結果

標準偏差を大きくしていくと、ぼかしの具合が強くなります。

カーネル=(15, 15)、標準偏差0

import cv2

img = cv2.imread('img/img.jpg')
dst = cv2.GaussianBlur(img, (15, 15), 0)
cv2.imwrite('img/img_gauss.jpg', dst)

結果

標準偏差に0を指定すると、カーネルサイズから標準偏差が自動計算されます。特に細かな調整が必要なければ、標準偏差は0にしてカネールサイズだけ変更すれば十分かと思います。

Python: OpenCVを使用して画像のサイズを変更する

OpenCVのインストール

pip install opencv-python

リサイズ

元の画像

画像の倍率を指定してリサイズ

import cv2

img = cv2.imread('img/img.jpg')
dst = cv2.resize(img, None, fx=0.5, fy=1.0)                                                                                     
cv2.imwrite('img/img_resized.jpg', dst)

実行結果

画像のサイズを指定してリサイズ

import cv2

img = cv2.imread('img/img.jpg')                                                                                                 
dst = cv2.resize(img, (300, 300))
cv2.imwrite('img/img_resized.jpg', dst)

実行結果

Python: Numpy入門

NumpyはPythonの数値計算用ライブラリです。Numpyを使用することで、行列計算を高速に行うことができます。

任意の1次元配列を作成

import numpy as	np
data = [1, 2, 3, 4]
arr = np.array(data)
print(arr)
#[1 2 3 4]
print(type(arr))
#<class 'numpy.ndarray'>
print(arr.shape)
#(4,)

任意の2次元配列を作成

import numpy as np
data = [[1, 2, 3, 4], [5, 6, 7, 8]]
arr = np.array(data)
print(arr)
#[[1 2 3 4]
# [5 6 7 8]]
print(arr.shape)
#(2, 4)

要素がすべて0の配列を作成

import numpy as np
arr = np.zeros(10)
print(arr)
#[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
print(arr.shape)
#(10,)
arr = np.zeros((3, 4))
print(arr)
#[[0. 0. 0. 0.]
# [0. 0. 0. 0.]
# [0. 0. 0. 0.]]
print(arr.shape)
#(3, 4)

要素がすべて1の配列を作成

import numpy as np
arr = np.ones(10)
print(arr)
#[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
print(arr.shape)
#(10,)
arr = np.ones((3, 4))
print(arr)
#[[1. 1. 1. 1.]
# [1. 1. 1. 1.]
# [1. 1. 1. 1.]]
print(arr.shape)
#(3, 4)

連続した数字の配列を作成

import numpy as	np
arr1 = np.arange(10)
print(arr1)
#[0 1 2 3 4 5 6 7 8 9]
print(arr1.shape)
#(10,)
arr2 = np.arange(5, 10)
print(arr2)
#[5 6 7 8 9]
print(arr2.shape)
#(5,)
arr3 = np.arange(2, 10, 2)
print(arr3)
#[2 4 6 8]
print(arr3.shape)
#(4,)

基本的な計算

import numpy as np
data1 = [[1, 2], [3, 4]]
arr1 = np.array(data1)
data2 = [[0, 1], [2, 3]]
arr2 = np.array(data2)

print(arr1)
#[[1 2]
# [3 4]]
print(arr1.shape)
#(2, 2)
print(arr2)
#[[0 1]
# [2 3]]
print(arr2.shape)
#(2, 2)

print(arr1 * arr2)
#[[ 0  2]
# [ 6 12]]
print(arr1 - arr2)
#[[1 1]
# [1 1]]
print(arr2 / arr1)
#[[0.         0.5       ]
# [0.66666667 0.75      ]]
print(arr1 ** 2)
#[[ 1  4]
# [ 9 16]]
print(arr1 > arr2)
#[[ True  True]
# [ True  True]]

Pythonの標準ライブラリdatetimeを使用して、日時の処理を行う

Pythonの標準ライブラリdatetimeを使用して、日時の処理を行います。

現在時刻の取得

import datetime
dt_now = datetime.datetime.now()
print(dt_now)
#2019-01-11 16:17:40.754123
print(type(dt_now))
#<class 'datetime.datetime'>
print(dt_now.year)
#2019
print(dt_now.month)
#1
print(dt_now.day)
#11
print(dt_now.hour)
#16
print(dt_now.minute)
#17
print(dt_now.second)
#40

任意のdatetimeオブジェクトを生成

dt = datetime.datetime(2019, 1, 1, 20, 00, 0)
print(dt_now)
#2019-01-11 16:23:39.727007
print(type(dt_now))
#<class 'datetime.datetime'>
print(dt_now.year)
#2019
print(dt_now.month)
#1
print(dt_now.day)
#11
print(dt_now.hour)
#16
print(dt_now.minute)
#23
print(dt_now.second)
#39

日時の時間差を取得

deltatimeオブジェクト同士を引き算することで、timedeltaオブジェクトが得られます。

dt_now = datetime.datetime.now()
print(dt_now)
#2019-01-11 16:29:47.220604
dt = datetime.datetime(2019, 1, 1, 20, 00, 0)
print(dt)
#2019-01-01 20:00:00
delta = dt_now - dt
print(delta)
#9 days, 20:29:47.220604
print(type(delta))
#<class 'datetime.timedelta'>
print(delta.days)
#9
print(delta.seconds)
#73787

datetimeオブジェクトを文字列に変換

dt_now = datetime.datetime.now()
dt_str = dt_now.strftime('%Y/%m/%d %H:%M')
print(dt_str)
#2019/01/11 16:34
print(type(dt_str))
#<class 'str'>

文字列をdatetimeオブジェクトに変換

dt = datetime.datetime.strptime('20180101', '%Y%m%d')
print(dt)
#2018-01-01 00:00:00
print(type(dt))
#<class 'datetime.datetime'>

Pythonでよく出るエラーメッセージとその解決方法

TypeError: can only concatenate str (not “int”) to str

error4.py
'1' + 1
実行
python error4.py 
Traceback (most recent call last):
  File "error4.py", line 1, in 
    '1' + 1
TypeError: can only concatenate str (not "int") to str

strとintは結合できません。

修正版
'1' + '1'

TypeError: unsupported operand type(s) for +: ‘int’ and ‘str’

error5.py
1 + '1'
実行
python error5.py 
Traceback (most recent call last):
  File "error5.py", line 1, in 
    1 + '1'
TypeError: unsupported operand type(s) for +: 'int' and 'str'

intとstrの結合はできません。

修正版
1 + 1

IndentationError: unexpected indent

error1.py
num = 1
    print(num)
実行
python error1.py 
  File "error1.py", line 2
    print(num)
    ^
IndentationError: unexpected indent

予期しないインデントが存在しています。

修正版
num = 1
print(num)

SyntaxError: EOL while scanning string literal

error2.py
print('aaa)
実行
python error2.py 
  File "error2.py", line 1
    print('aaa)
              ^
SyntaxError: EOL while scanning string literal

‘や”が閉じられていません。

修正版
print('aaa')

AttributeError: ‘int’ object has no attribute ‘append’

error3.py
mylist = 1
mylist.append(2)
実行
python error3.py 
Traceback (most recent call last):
  File "error3.py", line 2, in 
    mylist.append(2)
AttributeError: 'int' object has no attribute 'append'

「’int’ object」 = 「mylist」はappend属性を持っていません。

修正版
mylist = [1]
mylist.append(2)

TypeError: myfunc() missing 1 required positional argument: ‘arg2’

error6.py
def myfunc(arg1, arg2):
    print(arg1, arg2)

myfunc(1)
実行
python error6.py 
Traceback (most recent call last):
  File "error6.py", line 4, in 
    myfunc(1)
TypeError: myfunc() missing 1 required positional argument: 'arg2'

必要な引数’arg2’が指定されていません。

修正版
def myfunc(arg1, arg2):
    print(arg1, arg2)

myfunc(1, 2)

TypeError: myfunc() takes 2 positional arguments but 3 were given

error.py
def myfunc(arg1, arg2):
    print(arg1, arg2)

myfunc(1, 2, 3)
実行
python error7.py 
Traceback (most recent call last):
  File "error7.py", line 4, in 
    myfunc(1, 2, 3)
TypeError: myfunc() takes 2 positional arguments but 3 were given

引数2つの関数に3つの引数を与えています。

修正版
def myfunc(arg1, arg2):
    print(arg1, arg2)

myfunc(1, 2)

IndexError: list index out of range

error8.py
mylist = [1, 2, 3]
print(mylist[3])
実行
python error8.py 
Traceback (most recent call last):
  File "error8.py", line 2, in 
    print(mylist[3])
IndexError: list index out of range

指定したindexは存在しません。

修正版
mylist = [1, 2, 3]
print(mylist[2])

KeyError: ‘d’

error.py
mydict = {'a': 1, 'b': 2, 'c': 3}
print(mydict['d'])
実行
python error9.py 
Traceback (most recent call last):
  File "error9.py", line 2, in 
    print(mydict['d'])
KeyError: 'd'

指定したキーは存在しません。

修正版
mydict = {'a': 1, 'b': 2, 'c': 3}
print(mydict['c'])

Deep Learningを活用したゲノム変異解析ソフト「DeepVariant」を動かしてみる

「DeepVariant」はGoogle Brain と Verily Life Sciencesが開発したゲノム変異解析用のソフトウェアです。2016年には、ゲノム変異解析の精度を競う「PrecisionFDA Truth Challenge」で、現在デファクトスタンダードとなっている「GATK」などを抑えて「Highest SNP Performance賞」を受賞しました。「DeepVariant」は2017年12月4日にオープンソースとしてGitHub上にリリースされました。今回は、quick startに従って、Amazon EC2上で「DeepVariant」を実行してみます。

AWSコンソールから「サービス」->「EC2」->「インスタンス」->「インスタンスの作成」をクリックします。
Amazonマシンイメージには、「Ubuntu Server 18.04 LTS (HVM), SSD Volume Type – ami-0ac019f4fcb7cb7e6 (64 ビット x86)」、インスタンスタイプには「t2.micro」を選択しました。ストレージのサイズは8 GBで十分です。

まずは、dockerのインストールとイメージのダウンロードです。

BIN_VERSION="0.7.2"
MODEL_VERSION="0.7.2"

MODEL_NAME="DeepVariant-inception_v3-${MODEL_VERSION}+data-wgs_standard"
MODEL_HTTP_DIR="https://storage.googleapis.com/deepvariant/models/DeepVariant/${MODEL_VERSION}/${MODEL_NAME}"
DATA_HTTP_DIR="https://storage.googleapis.com/deepvariant/quickstart-testdata"

sudo apt -y update
sudo apt-get -y install docker.io
sudo docker pull gcr.io/deepvariant-docker/deepvariant:"${BIN_VERSION}"

続いて、Modelファイルをダウンロードしていきます。

mkdir -p ${MODEL_NAME}
wget -P ${MODEL_NAME} ${MODEL_HTTP_DIR}/model.ckpt.data-00000-of-00001
wget -P ${MODEL_NAME} ${MODEL_HTTP_DIR}/model.ckpt.index
wget -P ${MODEL_NAME} ${MODEL_HTTP_DIR}/model.ckpt.meta

[MODEL_NAME]フォルダの中に3つのファイルがダウンロードされたかと思います。これらは、TensorFlowのcheckpointフォーマットに当たります。

ls -1 "${MODEL_NAME}/"
#model.ckpt.data-00000-of-00001
#model.ckpt.index
#model.ckpt.meta

次に、サンプルデータをダウンロードします。

mkdir -p quickstart-testdata
wget -P quickstart-testdata "${DATA_HTTP_DIR}"/NA12878_S1.chr20.10_10p1mb.bam
wget -P quickstart-testdata "${DATA_HTTP_DIR}"/NA12878_S1.chr20.10_10p1mb.bam.bai
wget -P quickstart-testdata "${DATA_HTTP_DIR}"/test_nist.b37_chr20_100kbp_at_10mb.bed
wget -P quickstart-testdata "${DATA_HTTP_DIR}"/test_nist.b37_chr20_100kbp_at_10mb.vcf.gz
wget -P quickstart-testdata "${DATA_HTTP_DIR}"/test_nist.b37_chr20_100kbp_at_10mb.vcf.gz.tbi
wget -P quickstart-testdata "${DATA_HTTP_DIR}"/ucsc.hg19.chr20.unittest.fasta
wget -P quickstart-testdata "${DATA_HTTP_DIR}"/ucsc.hg19.chr20.unittest.fasta.fai
wget -P quickstart-testdata "${DATA_HTTP_DIR}"/ucsc.hg19.chr20.unittest.fasta.gz
wget -P quickstart-testdata "${DATA_HTTP_DIR}"/ucsc.hg19.chr20.unittest.fasta.gz.fai
wget -P quickstart-testdata "${DATA_HTTP_DIR}"/ucsc.hg19.chr20.unittest.fasta.gz.gzi

quickstart-testdataフォルダの中に10つのファイルがダウンロードされました。

ls -1 quickstart-testdata/
NA12878_S1.chr20.10_10p1mb.bam
NA12878_S1.chr20.10_10p1mb.bam.bai
test_nist.b37_chr20_100kbp_at_10mb.bed
test_nist.b37_chr20_100kbp_at_10mb.vcf.gz
test_nist.b37_chr20_100kbp_at_10mb.vcf.gz.tbi
ucsc.hg19.chr20.unittest.fasta
ucsc.hg19.chr20.unittest.fasta.fai
ucsc.hg19.chr20.unittest.fasta.gz
ucsc.hg19.chr20.unittest.fasta.gz.fai
ucsc.hg19.chr20.unittest.fasta.gz.gzi

NA12878_S1.chr20.10_10p1mb.bamは参照配列にマッピング済みのBAMファイル、ucsc.hg19.chr20.unittest.fastaは参照配列のFASTAファイル、test_nist.b37_chr20_100kbp_at_10mb.vcf.gzは評価用の正解データ、test_nist.b37_chr20_100kbp_at_10mb.bedは評価する領域を示したBEDファイルです。

それでは、実際に実行していきます。まず、結果出力用のフォルダ作成と環境変数の設定です。

OUTPUT_DIR=${HOME}/quickstart-output
mkdir -p "${OUTPUT_DIR}"
REF=${HOME}/quickstart-testdata/ucsc.hg19.chr20.unittest.fasta
BAM=${HOME}/quickstart-testdata/NA12878_S1.chr20.10_10p1mb.bam
MODEL="${HOME}/${MODEL_NAME}/model.ckpt"

STEP1: make_examples

このステップでは、pileup(ゲノム座標ごとに塩基をカウントしたもの)を作成します。

sudo docker run \
  -v ${HOME}:${HOME} \
  gcr.io/deepvariant-docker/deepvariant:"${BIN_VERSION}" \
  /opt/deepvariant/bin/make_examples \
  --mode calling   \
  --ref "${REF}"   \
  --reads "${BAM}" \
  --regions "chr20:10,000,000-10,010,000" \
  --examples "${OUTPUT_DIR}/examples.tfrecord.gz"

2つファイルが作成されました。

ls -1 quickstart-output/
examples.tfrecord.gz
examples.tfrecord.gz.run_info.pbtxt

STEP2: call_variants

このステップでは、実際に変異解析を実行します。

sudo docker run \
>   -v ${HOME}:${HOME} \
>   gcr.io/deepvariant-docker/deepvariant:"${BIN_VERSION}" \
>   /opt/deepvariant/bin/call_variants \
>  --outfile "${CALL_VARIANTS_OUTPUT}" \
>  --examples "${OUTPUT_DIR}/examples.tfrecord@${N_SHARDS}.gz" \
>  --checkpoint "${MODEL}"
I0110 02:48:42.534190 140314581722880 call_variants.py:292] Set KMP_BLOCKTIME to 0
Traceback (most recent call last):
  File "/tmp/Bazel.runfiles_hh_Ngv/runfiles/com_google_deepvariant/deepvariant/call_variants.py", line 411, in 
    tf.app.run()
  File "/root/.local/lib/python2.7/site-packages/tensorflow/python/platform/app.py", line 125, in run
    _sys.exit(main(argv))
  File "/tmp/Bazel.runfiles_hh_Ngv/runfiles/com_google_deepvariant/deepvariant/call_variants.py", line 401, in main
    use_tpu=FLAGS.use_tpu,
  File "/tmp/Bazel.runfiles_hh_Ngv/runfiles/com_google_deepvariant/deepvariant/call_variants.py", line 295, in call_variants
    example_format = tf_utils.get_format_from_examples_path(examples_filename)
  File "/tmp/Bazel.runfiles_hh_Ngv/runfiles/com_google_deepvariant/deepvariant/tf_utils.py", line 177, in get_format_from_examples_path
    one_example = get_one_example_from_examples_path(source)
  File "/tmp/Bazel.runfiles_hh_Ngv/runfiles/com_google_deepvariant/deepvariant/tf_utils.py", line 156, in get_one_example_from_examples_path
    'Cannot find matching files with the pattern "{}"'.format(source))
ValueError: Cannot find matching files with the pattern "/home/ubuntu/quickstart-output/examples.tfrecord@.gz"

何も考えずにコピペしたら、エラーになりました。–examples “${OUTPUT_DIR}/examples.tfrecord@${N_SHARDS}.gz”の部分が、並列で実行する場合になっていました。再度実行します。

sudo docker run \
  -v ${HOME}:${HOME} \
  gcr.io/deepvariant-docker/deepvariant:"${BIN_VERSION}" \
  /opt/deepvariant/bin/call_variants \
 --outfile "${CALL_VARIANTS_OUTPUT}" \
 --examples "${OUTPUT_DIR}/examples.tfrecord.gz" \
 --checkpoint "${MODEL}"

call_variants_output.tfrecord.gzが作成されました。

ls -1 quickstart-output/
call_variants_output.tfrecord.gz
examples.tfrecord.gz
examples.tfrecord.gz.run_info.pbtxt

STEP3: postprocess_variants

このステップでは変異検出結果をVCFファイルに変換します。

FINAL_OUTPUT_VCF="${OUTPUT_DIR}/output.vcf.gz"

sudo docker run \
  -v ${HOME}:${HOME} \
  gcr.io/deepvariant-docker/deepvariant:"${BIN_VERSION}" \
  /opt/deepvariant/bin/postprocess_variants \
  --ref "${REF}" \
  --infile "${CALL_VARIANTS_OUTPUT}" \
  --outfile "${FINAL_OUTPUT_VCF}"

output.vcf.gzが作成されました。

ls -1 quickstart-output/
call_variants_output.tfrecord.gz
examples.tfrecord.gz
examples.tfrecord.gz.run_info.pbtxt
output.vcf.gz

結果は以下のようなフォーマットになります。

##fileformat=VCFv4.2
##FILTER=
##FILTER=
##FILTER=
##INFO=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##contig=
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  NA12878
chr20   10000117        .       C       T       41.9    PASS    .       GT:GQ:DP:AD:VAF:PL      0/1:41:55:25,30:0.545455:41,0,49
chr20   10000211        .       C       T       46      PASS    .       GT:GQ:DP:AD:VAF:PL      0/1:46:59:30,29:0.491525:45,0,60
chr20   10000439        .       T       G       47.4    PASS    .       GT:GQ:DP:AD:VAF:PL      1/1:42:72:0,72:1:47,43,0
chr20   10000598        .       T       A       47.6    PASS    .       GT:GQ:DP:AD:VAF:PL      1/1:43:46:0,46:1:47,45,0
chr20   10000694        .       G       A       39.3    PASS    .       GT:GQ:DP:AD:VAF:PL      0/1:39:48:26,22:0.458333:39,0,63
chr20   10000758        .       T       A       46.5    PASS    .       GT:GQ:DP:AD:VAF:PL      1/1:43:56:0,56:1:46,46,0
chr20   10001019        .       T       G       4.6     PASS    .       GT:GQ:DP:AD:VAF:PL      0/1:5:44:31,13:0.295455:2,0,34
chr20   10001298        .       T       A       48.7    PASS    .       GT:GQ:DP:AD:VAF:PL      1/1:42:43:0,43:1:48,43,0
...

結果の評価

sudo docker pull pkrusche/hap.py
sudo docker run -it -v ${HOME}:${HOME} \
  pkrusche/hap.py /opt/hap.py/bin/hap.py \
  ${HOME}/quickstart-testdata/test_nist.b37_chr20_100kbp_at_10mb.vcf.gz \
  "${FINAL_OUTPUT_VCF}" \
  -f ${HOME}/quickstart-testdata/test_nist.b37_chr20_100kbp_at_10mb.bed \
  -r "${REF}" \
  -o "${OUTPUT_DIR}/happy.output" \
  --engine=vcfeval \
  -l chr20:10000000-10010000

良好な結果が得られました。

Benchmarking Summary:
  Type Filter  TRUTH.TOTAL  TRUTH.TP  TRUTH.FN  QUERY.TOTAL  QUERY.FP  QUERY.UNK  FP.gt  METRIC.Recall  METRIC.Precision  METRIC.Frac_NA  METRIC.F1_Score  TRUTH.TOTAL.TiTv_ratio  QUERY.TOTAL.TiTv_ratio  TRUTH.TOTAL.het_hom_ratio  QUERY.TOTAL.het_hom_ratio
 INDEL    ALL            4         4         0           13         0          9      0              1                 1        0.692308                1                     NaN                     NaN                   0.333333                   1.000000
 INDEL   PASS            4         4         0           13         0          9      0              1                 1        0.692308                1                     NaN                     NaN                   0.333333                   1.000000
   SNP    ALL           44        44         0           60         0         16      0              1                 1        0.266667                1                     1.2                1.307692                   0.333333                   0.363636
   SNP   PASS           44        44         0           60         0         16      0              1                 1        0.266667                1                     1.2                1.307692                   0.333333                   0.363636

AWS S3コマンドの使い方メモ

AWS Command Line Interfaceを使用してAWS S3を操作する方法をよく忘れるので、メモしておきます。

ファイルをダウンロード

aws s3 cp s3://[バケット名]/[ファイル名] .

ファイルをアップロード

aws s3 cp [ファイル名] s3://[バケット名]/

フォルダをダウンロード

aws s3 cp s3://[バケット名]/[フォルダ名]/ ./[フォルダ名] --recursive

フォルダをアップロード

aws s3 cp ./[フォルダ名] s3://[バケット名]/[フォルダ名]/ --recursive

Pythonで文字列処理をするときに、よく使うテクニック

文字列の結合

word1 = 'abc'
word2 = 'def'
word3 = word1 + '-' + word2
print(word3)
#abc-def
word1 = 'abc'
word2 = 'def'
word_list = [word1, word2]
print('-'.join(word_list))
#abc-def
word1 = 'abc'
word2 = 'def'
word3 = '{}-{}'.format(word1, word2)
print(word3)
#abc-def

文字列の切り出し

word = 'abcdefghijk'
print(word[3:])
#defghijk
print(word[:3])
#abc
print(word[1:3])
#bc
print(word[:-2])
#abcdefghi

文字列の置換

word = 'a-b-c-d-e'
new_word = word.replace('-', ',')
print(new_word)
#a,b,c,d,e

文字列の分割

word = 'a-b-c-d-e'
print(word.split('-'))
#['a', 'b', 'c', 'd', 'e']

文字列両端の空白削除

両端のスペース・タブ文字・改行を削除します

print(word.strip())
#'x'
word = ' x '
print(word.rstrip())
#' x'
print(word.lstrip())
#'x '