0

vcfファイルには最初の何行にも渡って#で書かれた箇所があり、読み込んだ時に#も一緒に出力してしまいます。その結果、counterでは{('3':987, '7':654, ・・・'#~~':1, '#~~':1)}といった感じに#のものまでカウントされてしまいます。この#を消す方法はあるのでしょうか?

また、カウントの順番を多い順ではなく、1,2,3と番号順にすることはできるのでしょうか?

import sys
import os
from collections import Counter

count = []
with open('test.vcf','r') as file:
    lines = file.read().split('\n')

    for line in lines:
        a = line.split('\t')
        CHR = a[0]
        count.append(CHR)

c = Counter(count)
print(c)
1
  • Bcftoolsを使う方法もあります(事前に VCF ファイルを bgzip で圧縮して index を作成しておく必要がありますが)。bcftools index --stats test.vcf.gz
    – user39889
    Commented 2019年5月15日 9:43

2 件の回答 2

1

先の回答で記載したとおりCounterクラスは辞書型のサブクラスであるため、これを利用して実装可能です。

具体的には、dict.keys()でキーを一覧し、それからdelを行うことで可能でしょう。

また、辞書は順序の概念がないためソートは不可能ですが、これもキーを任意の順序にソートしてアクセスすることで対応可能です。

c_keys = c.keys().sort()

for k in c_keys:
    if k[0] == '#':
        del c[k]
    else:
        print(k, c[k])
0

こんな感じで、基本のライブラリと pandas を使って実現できるでしょう。

import sys
import os
import pandas as pd

# ファイル読み取り&コメント削除(有効なデータの中に'#'は無いものと仮定)
lines = []
with open('test.vcf','r') as file:
    alltext = file.readlines()
    for line in alltext:
        str = line.strip().split('#')  # 行の前後の空白を削除し、'#'で分割
        if str[0]:                     # コメントでは無い部分だけリスト化
            lines.append(str[0])

# 各行のタブで区切った最初のデータをリスト化
count = []
for line in lines:
    a = line.split('\t')
    count.append(a[0])

# データの出現回数をカウントしてDataFrame化
s = pd.Series(count)
vc = pd.value_counts(s)
df = vc.rename_axis('CHR').to_frame('counts')

# データ文字列の最大長を取得
idx = list(df.index)
n = len(max(idx, key=len))

# 数値データを抽出し、右寄せしてソート
dfnum = df.query('CHR.str.isnumeric()', engine='python')
dfnum = dfnum.rename(index=lambda s: s.rjust(n))
dfnum = dfnum.sort_index()

# 文字列データを抽出し、右寄せしてソート
dfstr = df.query('not CHR.str.isnumeric()', engine='python')
dfstr = dfstr.rename(index=lambda s: s.rjust(n))
dfstr = dfstr.sort_index()

# 数値・文字列のDataFrameを連結
dfall = pd.concat([dfnum, dfstr])

# 結果表示
print(dfall)

検索したら、.vcf ファイル等を扱うライブラリがあるようですので、そちらを使うのもアリかもしれません。

cyvcf2 / pysam / PyVCF
cyvcf2: fast, flexible variant analysis with Python

他に古そうですがこんなまとめページも。
Python の基礎とバイオインフォマティクス Python 3.4


直接は関係無いですが、専門分野向けパッケージ等も。

Bioconda
biocondaを利用してNGS関連のソフトウェアを一括でインストールする
Biocondaを使ってみた

Python for Biologists

Biopython

2
  • この内容でPandasを使うのは牛刀という感じがします、VCFファイル専用のライブラリを使うのはありですが……
    – PicoSushi
    Commented 2019年5月17日 2:55
  • @PicoSushiさん、確かに import だけで二拍くらいかかっている感じです。けどその後の処理が素直に書ける感じなので良さそうかな、と思いまして。
    – kunif
    Commented 2019年5月17日 3:16

この質問に回答するには、ログインする必要があります。

求めていた回答ではありませんか? のタグが付いた他の質問を参照する。