鱼C论坛

 找回密码
 立即注册
查看: 1827|回复: 0

[分享] 请问这里错在哪里,为什么运行不了!FileNotFoundError: [Errno 2] No such file o...

[复制链接]
发表于 2021-1-28 01:19:03 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能^_^

您需要 登录 才可以下载或查看,没有账号?立即注册

x
import sys
import os
import subprocess

class Record(object):
    '''
    One line information in vcf file
    '''
    def __init__(self, line):
        info = line.split("\t")
        self.line = line
        self.CHROM =  info[0]
        self.POS = info[1]
        self.ID = info[2]
        self.REF = info[3]
        self.ALT = info[4]
        self.QUAL = info[5]
        self.FILTER = info[6]
        self.INFO = [{pair_lst[0]: pair_lst[1] if len(pair_lst)> 1 else ""} for pair_lst in [pair.split("=") for pair in info[7].split(";")]]
        self.FORMAT = info[8].split(":")
        self.sample_num = len(info) -7
        self.GT = []
        for i in range(2):
           GT_value = info[8 + i +1].split(":")
           GT_dict = {}
           for g in range(len(GT_value)):
               GT_dict[self.FORMAT[g]] = GT_value[g]
           self.GT.append(GT_dict)
      

class VCF(object):
    '''
    VCF class, read VCF, write VCF, get VCF information
    '''
    def __init__(self, uncompress_vcf):
        self.header = []
        self.reader = open(uncompress_vcf, 'r')
        self.line = self.reader.readline().strip()
        while self.line.startswith('#'):
            self.header.append(self.line)
            self.line = self.reader.readline().strip()
        self.record = Record(self.line)
    def __iter__(self):
        return self
    def __next__(self):
        self.line = self.reader.readline().strip()
        if self.line != "":
            self.record = Record(self.line)
            return self.record
        else:
            self.reader.close()
            raise StopIteration()
def reader_close(self):
        self.reader.close()   

gatk_result = "00007.vcf"
gatk = VCF(gatk_result)
gatk.header


record = gatk.record  #这里record存储的是该Record类的地址

record.CHROM
record.line
record.ID #其他的属性同理

CONTQ=28;DP=38;ECNT=1;GERMQ=76;MBQ=20,37;MFRL=171,229;MMQ=60,60;MPOS=26;NALOD=1.16;NLOD=3.91;
POPAF=6.00;RCNTS=0,0;ROQ=14;SEQQ=1;STRANDQ=11;TLOD=4.56
record.INFO

[{'CONTQ': '28'}, {'DP': '38'}, {'ECNT': '1'}, {'GERMQ': '76'}, {'MBQ': '20,37'}, {'MFRL': '171,229'}, {'MMQ': '60,60'}, {'MPOS': '26'}, {'NALOD': '1.16'}, {'NLOD': '3.91'}, {'POPAF': '6.00'}, {'RCNTS': '0,0'}, {'ROQ': '14'}, {'SEQQ': '1'}, {'STRANDQ': '11'}, {'TLOD': '4.56'}]

record.GT

[{'GT': '0/1', 'OBQRC': '41.80', 'SB': '12,9,1,1', 'DP': '23', 'OBF': '0.500', 'OBAM': 'false', 'OBP': '0.078', 'AD': '21,2', 'F2R1': '13,1', 'F1R2': '7,1', 'AF': '0.120', 'OBQ': '100.00', 'OBAMRC': 'true'}, {'GT': '0/0', 'OBQRC': '.', 'SB': '10,3,0,0', 'DP': '13', 'OBF': '.', 'OBAM': 'false', 'OBP': '.', 'AD': '13,0', 'F2R1': '6,0', 'F1R2': '7,0', 'AF': '0.065', 'OBQ': '.', 'OBAMRC': 'false'}]

for record in gatk:
    # compare GATK tumor AF to 0.05
    if float(record.GT[0]['AF']) > 0.05:
        print(record.line)

snv = "00007.vcf"

result = gatk.header
for record in gatk:
    if record.FILTER == "PASS" and float(record.GT[0]['AF']) > 0.05:
        result.append(record.line)

# write out result
with open(snv, 'w+') as snvf:
    for line in result:
        print(line, file = snvf)

record = next(gatk)
print(record.line)
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

小黑屋|手机版|Archiver|鱼C工作室 ( 粤ICP备18085999号-1 | 粤公网安备 44051102000585号)

GMT+8, 2024-5-17 12:02

Powered by Discuz! X3.4

© 2001-2023 Discuz! Team.

快速回复 返回顶部 返回列表