kuroの覚え書き

96の個人的覚え書き

homebrewにエラー再び

==> Upgrading ruby 
==> Downloading https://linuxbrew.bintray.com/bottles/ruby-2.6.1.x86_64_linux.bo
######################################################################## 100.0%
==> Pouring ruby-2.6.1.x86_64_linux.bottle.tar.gz
Error: The `brew link` step did not complete successfully
The formula built, but is not symlinked into /home/linuxbrew/.linuxbrew
Could not symlink bin/bundle
Target /home/linuxbrew/.linuxbrew/bin/bundle
is a symlink belonging to ruby. You can unlink it:
  brew unlink ruby

To force the link and overwrite all conflicting files:
  brew link --overwrite ruby

To list all files that would be deleted:
  brew link --overwrite --dry-run ruby

Possible conflicting files are:
/home/linuxbrew/.linuxbrew/bin/bundle -> /home/linuxbrew/.linuxbrew/Cellar/ruby/2.6.0_1/libexec/gembin/bundle
/home/linuxbrew/.linuxbrew/bin/bundler -> /home/linuxbrew/.linuxbrew/Cellar/ruby/2.6.0_1/libexec/gembin/bundler

んー。またrubyか。めんどくさいやっちゃなー。

あと、cufflinksがヘッドノードでは走るのに、計算ノードから起動しようとすると
Illegal inspiration
というエラーが出て起動しない。
これではtorqueでqsubできないのでめっちゃボトルネックができてしまう。

いろいろ試してみた。

$ brew remove cufflinks

で一旦アンインストールし、ヘッドノードで

$ brew install cufflinks

とインストールし直し→解決せず。
もう一度アンインストールし、計算ノードからインストール→これもだめ。
アンインストールしてから

$ brew cleanup

で不要なリンクファイルなどを一掃してからインストールし直し→解決。

今度からこまったときはbrew cleanupね。

CRISPRで編集がかかったゲノムシークエンスを読むプログラム

CRISPRで編集がかかったと思われるT0サンプルのゲノムをシークエンスすると、変異は2本の染色体にランダムに起こるため、基本的にはヘテロとなってシークエンスデータの波形が2つ分重なって検出されることになるのは先日に書いたとおり。
これをパズルのように分解して再構築して2本の染色体に分離し、それぞれの変異を同定しているわけだが、なんとか自動できないものかと。
というわけでプログラムを作ってみた。

ref.txt:入力としてreference配列のファイル
motif.txt:とりあえず2つのシークエンスが混ざっているのを記述したファイル
例えば
ref.txt
ATGGTGAGCATCAGCCTCTGGGTGG

ベースコールでは
ATGGTGAGCATGACTCTCTGTGG
となっていて、ATGGTGAGCまでは山がずれてなくて、以後2番手の山がかぶっている
ATGGTGAGCCACCGCCTGGGGGT
のように読み取れたとするなら
motif.txtに
ずれ始めから
A/CT/AG/CA/CC/GT/CCTC/GT/GGT/GGG/T
というふうに2つが重なっている部分は/で2つを併記する。

これを
crisgt.py

import sys
import argparse
import re

def get_args():
    parser = argparse.ArgumentParser()

    if sys.stdin.isatty():
        parser.add_argument("file", help="please set reference seq file", type=str)
        parser.add_argument("motif", help="please set search motif file", type=str)

    args = parser.parse_args()

    return(args)

def main():
    args = get_args()

    if hasattr(args, 'file'):
        file = args.file
    if hasattr(args, 'motif'):
        motif = args.motif

    with open(file) as list1:
        w = list1.read().replace('a','A').replace('c','C').replace('g','G').replace('t','T')
        with open(motif) as list2:
            k = list2.read().strip()
            k = k.replace('a','A').replace('c','C').replace('g','G').replace('t','T').replace('N', '[ACGT]').replace('n', '[ACGT]').replace('A/T','[AT]').replace('A/C','[AC]').replace('A/G','[AG]').replace('C/G','[CG]').replace('C/T','[CT]').replace('G/T','[GT]')
            m = re.finditer(k,w)
            for match in m:
               print(match.group(), match.start(), match.end())

if __name__ == '__main__':
    main()
$ python3 crisgt.py ref.txt motif.txt
ATCAGCCTCTGGGT 9 23

というふうに処理すると
あれれ。2つ目の位置が記録されないぞ?

qiita.com

お、これか。pythonのmatchは重複を認めない設定になっているのか。こりゃ以前に作ったプログラムも見直さないといかんな。

            m = re.finditer(k,w)
import regex
            m = regex.finditer(k,w,overlapped=True)

これで解決。

$ python3 crisgtype.py ref.txt motif.txt 
ATCAGCCTCTGGGT 9 23
CAGCCTCTGGGTGG 11 25

これにより、9(10塩基め)から正常な配列に対し、2塩基抜けて12塩基目がずれて来ていることがわかる。

あとは重なった波形から自動的に2つのピークを抽出してきたりしてくれるとすごくいいが、そこはちょっとむずかしいかも。

import sys
import argparse
import re
import regex

def get_args():
    parser = argparse.ArgumentParser()

    if sys.stdin.isatty():
        parser.add_argument("file", help="please set reference seq file", type=str)
        parser.add_argument("seq1", help="please set seq1 file", type=str)
        parser.add_argument("seq2", help="please set seq2 file", type=str)

    args = parser.parse_args()

    return(args)

def main():
    args = get_args()

    if hasattr(args, 'file'):
        file = args.file
    if hasattr(args, 'seq1'):
        seq1 = args.seq1
    if hasattr(args, 'seq2'):
        seq2 = args.seq2

    with open(file) as list1:
        w = list1.read().replace('a','A').replace('c','C').replace('g','G').replace('t','T')
        with open(seq1) as s1:
            sq1 = s1.read().strip()
        with open(seq2) as s2:
            sq2 = s2.read().strip()
        s = "".join('['+i+j+']' for i, j in zip(sq1, sq2))
        m = regex.finditer(s,w,overlapped=True)

        for match in m:
            print(match.group(), match.start(), match.end())

if __name__ == '__main__':
    main()

お、ずいぶんスッキリしてきたぞ。
山が2つ重なりだすところからのベースコールデータと、2番手の山を読み取った配列をそれぞれseq1.txt、seq2.txtとしておいてこれを

$ python3 crisgt.py ref.txt seq1.txt seq2.txt
ATCAGCCTCTGGGT 9 23
CAGCCTCTGGGTGG 11 25

こりゃもう実用レベルじゃね?

CRISPRでknockoutを作ったときにgenotypingをsanger sequenceすることで行う

CRISPRでINDELを誘発したゲノム配列が、実際どういうふうに編集されたかを確認するのにわざわざNGSをつかうのはちょっと大げさなので普通のsangerシークエンサーでシークエンスを読んで確認をしたい。しかし、変異は普通ヘテロに入るので、変異が入った部分からシークエンスの波形がずれて重なり合ってくる。
なので、2つの山を分離して、それぞれのシークエンスを読みほぐす必要がある。まあ、波形を表示するアプリケーションで開いて山をたどりながらエディタでACGTを書いていけばいいだけなんだけど、これが結構面倒くさい。

きっと同じようなことを考えてこれを自動でやっちゃうソフトを作っている人がいるはずだ、と思って検索してみたところ、3つ見つかった。

https://tide.nki.nl

Synthego

CRISP-ID: Detecting CRISPR mediated indels by Sanger sequencing


しかしどれもweb tool。まあいまどきだいたいそんな感じなんだろうけど、ネットに研究データを安直に流すのはどうよ、ってことでローカルでできるアプリがないかと探してみたところ、2つ目のICEについてはオープンソースで公開しているから、ローカルにインストールもできちゃうようだ。それもpythonっぽいので好都合。

ということで

$ git clone git@github.com/synthego-open/ice.git
fatal: repository 'git@github.com/synthego-open/ice.git' does not exist

おや?

$ git clone https://github.com/synthego-open/ice
Cloning into 'ice'...
remote: Enumerating objects: 336, done.
remote: Total 336 (delta 0), reused 0 (delta 0), pack-reused 336
Receiving objects: 100% (336/336), 2.02 MiB | 1.68 MiB/s, done.
Resolving deltas: 100% (203/203), done.

こっちでOKね。

$ cd ice
$ pip3 install -r requirements.txt
Collecting pytest==3.4.2 (from -r requirements.txt (line 1))
  Downloading https://files.pythonhosted.org/packages/f1/5c/411ceafef3b5e5486d16f174db18dc26f49e7704dbf59ef488e95db47339/pytest-3.4.2-py2.py3-none-any.whl (189kB)
    100% |████████████████████████████████| 194kB 4.2MB/s 
Requirement already satisfied: biopython>=1.70 in /Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages (from -r requirements.txt (line 2)) (1.70)
Collecting coverage>=4.4.1 (from -r requirements.txt (line 3))
  Downloading https://files.pythonhosted.org/packages/fb/af/ce7b0fe063ee0142786ee53ad6197979491ce0785567b6d8be751d2069e8/coverage-4.5.2.tar.gz (384kB)
    100% |████████████████████████████████| 389kB 4.2MB/s 
Collecting pytest-cov==2.5.1 (from -r requirements.txt (line 4))
  Downloading https://files.pythonhosted.org/packages/30/7d/7f6a78ae44a1248ee28cc777586c18b28a1df903470e5d34a6e25712b8aa/pytest_cov-2.5.1-py2.py3-none-any.whl
Collecting matplotlib>=2.1.0 (from -r requirements.txt (line 5))
  Downloading https://files.pythonhosted.org/packages/28/6c/addb3560777f454b1d56f0020f89e901eaf68a62593d4795e38ddf24bbd6/matplotlib-3.0.2-cp36-cp36m-macosx_10_6_intel.macosx_10_9_intel.macosx_10_9_x86_64.macosx_10_10_intel.macosx_10_10_x86_64.whl (14.1MB)
    100% |████████████████████████████████| 14.1MB 892kB/s 
Collecting numpy>=1.13.3 (from -r requirements.txt (line 6))
  Downloading https://files.pythonhosted.org/packages/88/b8/569d9c702685b595812fbfd9ee04f240653b7a15feec43cc98be3b34e5f5/numpy-1.16.1-cp36-cp36m-macosx_10_6_intel.macosx_10_9_intel.macosx_10_9_x86_64.macosx_10_10_intel.macosx_10_10_x86_64.whl (13.9MB)
    100% |████████████████████████████████| 13.9MB 1.5MB/s 
Collecting pandas>=0.20.3 (from -r requirements.txt (line 7))
  Downloading https://files.pythonhosted.org/packages/99/12/bf4c58eea94cea4f91ff931f284146337814fb8546e6eb0b52584446fd52/pandas-0.24.1-cp36-cp36m-macosx_10_6_intel.macosx_10_9_intel.macosx_10_9_x86_64.macosx_10_10_intel.macosx_10_10_x86_64.whl (16.3MB)
    100% |████████████████████████████████| 16.3MB 732kB/s 
Requirement already satisfied: scipy>=0.19.1 in /Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages (from -r requirements.txt (line 8)) (0.19.1)
Collecting sklearn>=0.0 (from -r requirements.txt (line 9))
  Downloading https://files.pythonhosted.org/packages/1e/7a/dbb3be0ce9bd5c8b7e3d87328e79063f8b263b2b1bfa4774cb1147bfcd3f/sklearn-0.0.tar.gz
Collecting xlrd>=1.1.0 (from -r requirements.txt (line 10))
  Downloading https://files.pythonhosted.org/packages/b0/16/63576a1a001752e34bf8ea62e367997530dc553b689356b9879339cf45a4/xlrd-1.2.0-py2.py3-none-any.whl (103kB)
    100% |████████████████████████████████| 112kB 10.0MB/s 
Collecting xlsxwriter>=1.0.2 (from -r requirements.txt (line 11))
  Downloading https://files.pythonhosted.org/packages/3d/1b/4caecd4efde1d41ba3bef1a81027032a7a6dff7d5112e1731f232c0addb9/XlsxWriter-1.1.2-py2.py3-none-any.whl (142kB)
    100% |████████████████████████████████| 143kB 370kB/s 
Requirement already satisfied: attrs>=17.2.0 in /Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages (from pytest==3.4.2->-r requirements.txt (line 1)) (18.2.0)
Requirement already satisfied: py>=1.5.0 in /Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages (from pytest==3.4.2->-r requirements.txt (line 1)) (1.7.0)
Collecting pluggy<0.7,>=0.5 (from pytest==3.4.2->-r requirements.txt (line 1))
  Downloading https://files.pythonhosted.org/packages/ba/65/ded3bc40bbf8d887f262f150fbe1ae6637765b5c9534bd55690ed2c0b0f7/pluggy-0.6.0-py3-none-any.whl
Requirement already satisfied: setuptools in /Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages (from pytest==3.4.2->-r requirements.txt (line 1)) (28.8.0)
Requirement already satisfied: six>=1.10.0 in /Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages (from pytest==3.4.2->-r requirements.txt (line 1)) (1.10.0)
Collecting kiwisolver>=1.0.1 (from matplotlib>=2.1.0->-r requirements.txt (line 5))
  Downloading https://files.pythonhosted.org/packages/fb/96/619db9bf08f652790fa9f3c3884a67dc43da4bdaa185a5aa2117eb4651e1/kiwisolver-1.0.1-cp36-cp36m-macosx_10_6_intel.macosx_10_9_intel.macosx_10_9_x86_64.macosx_10_10_intel.macosx_10_10_x86_64.whl (108kB)
    100% |████████████████████████████████| 112kB 2.4MB/s 
Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages (from matplotlib>=2.1.0->-r requirements.txt (line 5)) (2.2.0)
Requirement already satisfied: cycler>=0.10 in /Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages (from matplotlib>=2.1.0->-r requirements.txt (line 5)) (0.10.0)
Requirement already satisfied: python-dateutil>=2.1 in /Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages (from matplotlib>=2.1.0->-r requirements.txt (line 5)) (2.6.1)
Requirement already satisfied: pytz>=2011k in /Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages (from pandas>=0.20.3->-r requirements.txt (line 7)) (2017.2)
Collecting scikit-learn (from sklearn>=0.0->-r requirements.txt (line 9))
  Downloading https://files.pythonhosted.org/packages/cb/5f/dfa0a118b8a503e45cd2cf48acb9cf1de8deaf06a3cef1b1c19bd5cbbc45/scikit_learn-0.20.2-cp36-cp36m-macosx_10_6_intel.macosx_10_9_intel.macosx_10_9_x86_64.macosx_10_10_intel.macosx_10_10_x86_64.whl (7.9MB)
    100% |████████████████████████████████| 7.9MB 2.5MB/s 
Installing collected packages: pluggy, pytest, coverage, pytest-cov, kiwisolver, numpy, matplotlib, pandas, scikit-learn, sklearn, xlrd, xlsxwriter
  Found existing installation: pluggy 0.8.1
    Uninstalling pluggy-0.8.1:
      Successfully uninstalled pluggy-0.8.1
  Found existing installation: pytest 4.2.0
    Uninstalling pytest-4.2.0:
      Successfully uninstalled pytest-4.2.0
  Running setup.py install for coverage ... done
  Found existing installation: numpy 1.13.0
    Uninstalling numpy-1.13.0:
      Successfully uninstalled numpy-1.13.0
  Found existing installation: matplotlib 2.0.2
    Uninstalling matplotlib-2.0.2:
      Successfully uninstalled matplotlib-2.0.2
  Running setup.py install for sklearn ... done
Successfully installed coverage-4.5.2 kiwisolver-1.0.1 matplotlib-3.0.2 numpy-1.16.1 pandas-0.24.1 pluggy-0.6.0 pytest-3.4.2 pytest-cov-2.5.1 scikit-learn-0.20.2 sklearn-0.0 xlrd-1.2.0 xlsxwriter-1.1.2

でけた。

使い方は

Running a single ICE analysis:

./ice_analysis_single.py \
	--control ./ice/tests/test_data/good_example_control.ab1  \
	--edited ./ice/tests/test_data/good_example_edited.ab1 \
	--target AACCAGTTGCAGGCGCCCCA \
	--out results/testing \
	--verbose
Running a batch analysis:

./ice_analysis_batch.py \
	--in ./ice/tests/test_data/batch_example.xlsx \
	--out ./results/ \
	--data ./ice/tests/test_data/
	--verbose

だそうな。
やってみる

$ ./ice_analysis_single.py --control ./ice/tests/test_data/good_example_control.ab1 --edited ./ice/tests/test_data/good_example_edited.ab1 --target AACCAGTTGCAGGCGCCCCA --out results/testing --verbose
Synthego ICE (https://synthego.com)
Version: 1.1.1
Base dir: /Users/kkuro/local/bin/ice/results
analyzing 462 number of edit proposals
Shape of coefficient matrix: (462, 228)
------------------------------------------------------
Inference Sequence length: 57
Inference sequence 
Output vector : 4x57 (228)

NNLS input shapes
--------------------------------------------------------
A (228, 462)
b (228,)
R_SQUARED 0.9823262526685657
discord (aln window): 0.14 after cutsite: 0.60

結果は

control                             NNN-NN-NGGGGCATCCTGTGTTCTACCTGGCACCTGTCCCCATAGAAAT
edited                              NNNANNANNGGGCATCCTGTGTTCTACCTGGCACCTGTCCCCATAGAAAT

control                             GAGCGTGAGTGCCCGGGATCTGCTGCGGGGCTGTGCTGGGCTCTTTCTCA
edited                              GAGCGTGAGTGCCCGGGATCTGCTGCGGGGCTGTGCTGGGCTCTTTCTCA

control                             GCCTGGCCCGAAGTTTCCAGATCTGATTGAGCGAGAGAGCAGCAGGACCT
edited                              GCCTGGCCCGAAGTTTCCAGATCTGATTGAGCGAGAGAGCAGCAGGACCT

control                             GCCCCTCTGCTGGGCTCTTACCTTCGCGGCACTCGCCACTGCCCAGCAGC
edited                              GCCCCTCTGCTGGGCTCTTACCTTCGCGGCACTCGCCACTGCCCAGCAGC

control                             AGGTGAGGCCCAACACAACCAGTTGCAGGCGCCCCATGGTGAGCATCAGC
edited                              AGGTGAGGCCCAACACAACCAGTTGCAGGCGCCCC-TGGGGAGAATCACC

control                             CTCTGGGTGGCCCTCCCTCTGGGCCTCGGGTATTTATGGAGCTGGATCCA
edited                              CCCTGGGGGGCCCCCCCCCTGGGGCCCGGGGATTTTTGGGGAGGGAGCCC

control                             AGGTCACATGCTTGTTCATGAGCTCTC-AGGCA-
edited                              AGGGGCCATGTTTTTTTTTAAACCCCCNAGAAAA

こんな感じだね。ターゲットのPAM配列TGGの直前のAがdeleteされていることが検出されているな。
シークエンスデータは

f:id:k-kuro:20190204191845p:plain
シークエンス波形データ
こんな感じだった。
コントロールも波形データでないといかんのかなあ→fasta.txtだとだめだった。むぅ。

塩基配列(文字列)のどこにリストに挙げた配列が位置しているかを示すプログラム

モチーフの抽出ができたらそれが一体どこにあるのか列記したい。

import re

text = """GACTTTAGATGGCTTCTTCCTTTATAACCAATTGATATTGCATACTCTGATGAGATTTAT 
AATTAAAGAAGCAGAAACAAAAACAAGTAAAACAGAAACAATCAACACAGAGAAACCACC 
CCGAGAATATCTCCATTGGTTGGTGACTGATATCCCTGCTACAACTGGAACAACCTTTGG"""

list1 = text.split()

motif = """AAACA AAACAA AACAA AGAAA AGAAAC GAAAC"""

list2= motif.split()

for i, w in enumelate(list1):
    for k in list2:
        m = re.finditer(k,w)
        for match in m:
            print(i, match)

実行すると

$ python3 motif1.py
1 <_sre.SRE_Match object; span=(14, 19), match='AAACA'>
1 <_sre.SRE_Match object; span=(20, 25), match='AAACA'>
1 <_sre.SRE_Match object; span=(29, 34), match='AAACA'>
1 <_sre.SRE_Match object; span=(35, 40), match='AAACA'>
1 <_sre.SRE_Match object; span=(14, 20), match='AAACAA'>
1 <_sre.SRE_Match object; span=(20, 26), match='AAACAA'>
1 <_sre.SRE_Match object; span=(35, 41), match='AAACAA'>
1 <_sre.SRE_Match object; span=(15, 20), match='AACAA'>
1 <_sre.SRE_Match object; span=(21, 26), match='AACAA'>
1 <_sre.SRE_Match object; span=(36, 41), match='AACAA'>
1 <_sre.SRE_Match object; span=(12, 17), match='AGAAA'>
1 <_sre.SRE_Match object; span=(33, 38), match='AGAAA'>
1 <_sre.SRE_Match object; span=(50, 55), match='AGAAA'>
1 <_sre.SRE_Match object; span=(12, 18), match='AGAAAC'>
1 <_sre.SRE_Match object; span=(33, 39), match='AGAAAC'>
1 <_sre.SRE_Match object; span=(50, 56), match='AGAAAC'>
1 <_sre.SRE_Match object; span=(13, 18), match='GAAAC'>
1 <_sre.SRE_Match object; span=(34, 39), match='GAAAC'>
1 <_sre.SRE_Match object; span=(51, 56), match='GAAAC'>
2 <_sre.SRE_Match object; span=(48, 53), match='AACAA'>

こんな感じ。
span=()のところに存在位置が表示されている
このあとの処理がやりやすいようにもうちょっと整形しておきたいが。

import re

text = """GACTTTAGATGGCTTCTTCCTTTATAACCAATTGATATTGCATACTCTGATGAGATTTAT
AATTAAAGAAGCAGAAACAAAAACAAGTAAAACAGAAACAATCAACACAGAGAAACCACC
CCGAGAATATCTCCATTGGTTGGTGACTGATATCCCTGCTACAACTGGAACAACCTTTGG"""

list1 = text.split()

motif = """AAACA AAACAA AACAA AGAAA AGAAAC GAAAC"""

list2= motif.split()

for i, w in enumerate(list1):
    for k in list2:
        m = re.finditer(k,w)
        for match in m:
            print(i, match.group(), match.start())

これだ。

$ python3 motif2.py
1 AAACA 14
1 AAACA 20
1 AAACA 29
1 AAACA 35
1 AAACAA 14
1 AAACAA 20
1 AAACAA 35
1 AACAA 15
1 AACAA 21
1 AACAA 36
1 AGAAA 12
1 AGAAA 33
1 AGAAA 50
1 AGAAAC 12
1 AGAAAC 33
1 AGAAAC 50
1 GAAAC 13
1 GAAAC 34
1 GAAAC 51
2 AACAA 48

なお、match.span()とすると(start, end)のタプルが返る。end()だとendだけ。

検索したい配列とモチーフをそれぞれファイルで用意しておいて(各1行に1配列)スクリプトに食わせる

import sys
import argparse
import re

def get_args():
    parser = argparse.ArgumentParser()

    if sys.stdin.isatty():
        parser.add_argument("file", help="please set seq file", type=str)
        parser.add_argument("motif", help="please set motif list", type=str)

    args = parser.parse_args()

    return(args)

def main():
    args = get_args()

    if hasattr(args, 'file'):
        file = args.file
    if hasattr(args, 'motif'):
        motif = args.motif

    with open(file) as list1:
        for i, w in enumerate(list1):
            with open(motif) as list2:
                for k in list2:
                    m = re.finditer(k.strip(),w)
                    for match in m:
                        print(i, match.group(), match.start())

if __name__ == '__main__':
    main()

第一引数に配列ファイル、第二引数にモチーフファイルをつけて実行。出力はそのままだと標準出力。


2/11追記
reのfinditerでは重複したモチーフは飛ばされてしまうことが判明したので修正

import sys
import argparse
import regex

def get_args():
    parser = argparse.ArgumentParser()

    if sys.stdin.isatty():
        parser.add_argument("file", help="please set seq file", type=str)
        parser.add_argument("motif", help="please set motif list", type=str)

    args = parser.parse_args()

    return(args)

def main():
    args = get_args()

    if hasattr(args, 'file'):
        file = args.file
    if hasattr(args, 'motif'):
        motif = args.motif

    with open(file) as list1:
        for i, w in enumerate(list1):
            with open(motif) as list2:
                for k in list2:
                    m = regex.finditer(k.strip(),w,overlapped=True)
                    for match in m:
                        print(i, match.group(), match.start(), match.end())

if __name__ == '__main__':
    main()

これでOKのはず。

画像の白飛びをなんとかしたい

カメラの露出制御はマニュアルで行えるらしい。

Atelier Orchard: ホワイトバランス設定してウェブカメラ録画


なるほど。
というわけでmotionの画面を見ながら明るさやコントラストの調節を手動でできるようにアプリを作成。

f:id:k-kuro:20190128191008p:plain

こんな感じにコントロールできるようにアプリを作り込んだ。
業務で使う分にはここまで完成度を上げる必要もないのだが、これなら誰でも使えるだろう。

なお、v4lの設定項目は接続するカメラごとに最大値、最小値がちがうし、そもそも設定できる項目も違っているらしいい。なのでカメラを変えたら、プログラムも書き換えが必要になってくる。ちょっと面倒だな

そもそもタイムラプスの撮影間隔とかcrontabで制御している時点で誰でも設定して使える仕様ではないし、万人が使える装置とするにはまだまだ完成形とは言えないね。それとRTCをつけないとネットに繋がっていない環境で日時をちゃんと設定できないという問題もあるな。
って、売るわけじゃないからそこまでやらなくてもいいだろうけど。

avconv (ffmpeg)でmp4ビデオを作成してwebに埋め込む

いやーはまったはまった。
ラズパイ3で作ったwebカメラアプリを3+に入れてやったところ、ほんの数ヶ月の間にソフトウェアの構成がかなり変わったらしく、いろいろ言うことを聞かない。
最後まで手こずったのはタイムラプスで撮った写真をパラパラ漫画的にビデオに作り、できたビデオをwebに貼り付けてプレビューできるようにしていたところ。3+で作成されるビデオがMacsafariiOSで再生できない。
これではiPhoneで作業が完結できないじゃないか。

avconvという名前で使っていたffmpegの派生ソフトが元の鞘に収まってffmpegに戻った際に、どうもいろいろな初期設定が変わった(というかffmpegの初期設定に戻った)ようなのだ。

Why won't video from ffmpeg show in QuickTime, iMovie or quick preview? - Ask Different
ここの書き込みが決定打となってなんとか表示ができるようになった。

つまり

In short, you (often) need to include the argument -pix_fmt yuv420p when using ffmpeg to generate H.264 content for Apple software/devices, and a bunch of other decoders that don't handle yuv444p.

ここ。

-pix_fmt yuv420p

をoutputのオプションに入れてやったところ、無事に再生できる動画が生成されるようになった。

連番で同じ処理をさせるためのリストを作る

ファイルを連番を付けて作成していてそれに同じ操作を全部やるとき、excelで編集してコマンドリストを作っていたが、shellでやるほうが簡単なので、その覚え書き

seq -w 0 42 | xargs -i echo "python3 seq.py TAIR10_cdna_20101214_updated.fa.split/TAIR10_cdna_20101214_updated.part_0{}.fasta > TAIR10_{}.txt" > com.txt

たとえばこのように42分割したfastaファイルにそれぞれある処理を実行する。

塩基配列(文字列)からコンセンサス配列を抽出するプログラム

要するに文章からよく出てくる単語をピックアップしてカウントし、リストを作れれば良い。
辞書型を使って単語を数えるプログラムはpythonのプログラム例としてよく上がっているが、単語の区切りが明確でない遺伝子配列のような文字列から指定文字数の連続した文字列が指定反復以上出現していたらピックアップしてカウントするというプログラムは落ちていなかったので、作ってみた。多分日本語みたいに単語ごとに明確に区切りを入れずに書く文章から抽出したりするのにも使えると思う。その場合いわゆるキーワード抽出というのとは違い、文字を単純に検索しているだけなので、日本語としておかしな部分が切り取られる可能性も大いにあるが。

text = """ここに検索したい文字列(遺伝子配列など)を入れる。複数の文字列から共通する文字列を検索したいときはスペースで区切って複数並べておく"""
low = 4 #ここで指定した文字数以上の連続したコンセンサスを探す
time = 3 #最低反復数

list = text.split()
counter = {}
for w in list:
    l = len(w)
    l1 = l - low + 1
    for i in range(0,l1):
        for n in range(0,l1-i):
            ws = w[i:i+low+n]
            if ws in counter:
                counter[ws] += 1
            else:
                counter[ws] = 1
for k,v in sorted(counter.items()):
    if v >= time:
        p = text.find(k) #最初に出現する位置
        print(k,v,p)

ウェブインターフェースとか使いやすくする事はできるが、面倒臭いので生スクリプトで済ます。

homebrewにエラー

いつの頃からか

$ brew update
fatal: unable to access 'https://github.com/brewsci/homebrew-bio/': error setting certificate verify locations:
  CAfile: /home/linuxbrew/.linuxbrew/etc/openssl/cert.pem
  CApath: /home/linuxbrew/.linuxbrew/etc/openssl/certs
fatal: unable to access 'https://github.com/Linuxbrew/homebrew-core/': error setting certificate verify locations:
  CAfile: /home/linuxbrew/.linuxbrew/etc/openssl/cert.pem
  CApath: /home/linuxbrew/.linuxbrew/etc/openssl/certs
fatal: unable to access 'https://github.com/Linuxbrew/brew/': error setting certificate verify locations:
  CAfile: /home/linuxbrew/.linuxbrew/etc/openssl/cert.pem
  CApath: /home/linuxbrew/.linuxbrew/etc/openssl/certs
fatal: unable to access 'https://github.com/brewsci/homebrew-science/': error setting certificate verify locations:
  CAfile: /home/linuxbrew/.linuxbrew/etc/openssl/cert.pem
  CApath: /home/linuxbrew/.linuxbrew/etc/openssl/certs
fatal: unable to access 'https://github.com/Linuxbrew/homebrew-extra/': error setting certificate verify locations:
  CAfile: /home/linuxbrew/.linuxbrew/etc/openssl/cert.pem
  CApath: /home/linuxbrew/.linuxbrew/etc/openssl/certs
fatal: unable to access 'https://github.com/Linuxbrew/homebrew-xorg/': error setting certificate verify locations:
  CAfile: /home/linuxbrew/.linuxbrew/etc/openssl/cert.pem
  CApath: /home/linuxbrew/.linuxbrew/etc/openssl/certs
Error: Fetching /home/linuxbrew/.linuxbrew/Homebrew/Library/Taps/brewsci/homebrew-bio failed!
Fetching /home/linuxbrew/.linuxbrew/Homebrew/Library/Taps/homebrew/homebrew-core failed!
Fetching /home/linuxbrew/.linuxbrew/Homebrew failed!
Fetching /home/linuxbrew/.linuxbrew/Homebrew/Library/Taps/brewsci/homebrew-science failed!
Fetching /home/linuxbrew/.linuxbrew/Homebrew/Library/Taps/linuxbrew/homebrew-extra failed!
Fetching /home/linuxbrew/.linuxbrew/Homebrew/Library/Taps/linuxbrew/homebrew-xorg failed!

というようなエラーが出るようになった。brewで入れたアプリの挙動も若干おかしいような感じもあった。

続きを読む

motion.confの改変内容メモ

############################################################
# Daemon
############################################################

# Start in daemon (background) mode and release terminal (default: off)
daemon on

###########################################################
# Capture device options
############################################################

# Image width (pixels). Valid range: Camera dependent, default: 352
width 640

# Image height (pixels). Valid range: Camera dependent, default: 288
height 480

############################################################
# Image File Output
############################################################

# Output 'normal' pictures when motion is detected (default: on)
# Valid values: on, off, first, best, center
# When set to 'first', only the first picture of an event is saved.
# Picture with most motion of an event is saved when set to 'best'.
# Picture with motion nearest center of picture is saved when set to 'center'.
# Can be used as preview shot for the corresponding movie.
output_pictures off

############################################################
# Live Stream Server
############################################################

# Restrict stream connections to localhost only (default: on)
stream_localhost off