使用linux或python从文件中提取特定的列和字符串

问题描述 投票：0回答：2

处理12 Gb文件时遇到问题。谈到linux，我是新手。我希望有人能帮助我，任何建议都表示赞赏。

我有一个名为phase_3.vcf的文件，如下所示：

##INFO=<ID=EAS_AF,Number=A,Type=Float,Description="Allele frequency in the EAS populations">                            
##INFO=<ID=EUR_AF,Number=A,Type=Float,Description="Allele frequency in the EUR populations">                            
##INFO=<ID=AFR_AF,Number=A,Type=Float,Description="Allele frequency in the AFR populations">                            
##INFO=<ID=AMR_AF,Number=A,Type=Float,Description="Allele frequency in the AMR populations">                            
##INFO=<ID=SAS_AF,Number=A,Type=Float,Description="Allele frequency in the SAS populations">                
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO
1   10177   rs367896724 A   AC  .   .    dbSNP_150;E_Freq;E_1000G;EAS_AF=0.3363;SAS_AF=0.4949;AFR_AF=0.4909
1   10235   rs540431307 T   TA  .   .   dbSNP_150;E_Freq;E_1000G;EAS_AF=0.0000;AMR_AF=0.0014;
1   10352   rs555500075 T   TA  .   .   dbSNP_150;EAS_AF=0.4306;EUR_AF=0.4264;SAS_AF=0.4192;
1   10505   rs548419688 A   T   .   .   TSA=SNV;E_Freq;EAS_AF=0.0000;AMR_AF=0.0000;AFR_AF=0.0008
1   10506   rs568405545 C   G   .   .   dbSNP_150;TSA=SNV;MA=G;MAF=0.000199681;EAS_AF=0.0000;

我想保留前5列和字符串“EAS_AF =”，数字跟随它们。

为简单起见，名为result.txt的结果的预期形式应如下所示：

#CHROM  POS ID  REF ALT INFO
1   10177   rs367896724 A   AC  EAS_AF=0.3363
1   10235   rs540431307 T   TA  EAS_AF=0.0000
1   10352   rs555500075 T   TA  EAS_AF=0.4306
1   10505   rs548419688 A   T   EAS_AF=0.0000
1   10511   rs534229142 G   A   EAS_AF=0.0000
1   10539   rs537182016 C   A   EAS_AF=0.0000
python linux string pattern-matching bioinformatics
2个回答
2
投票

您的文件采用非常常见的生物信息学格式（vcf）。因此，生物信息学中有一些工具专门用于解决您的任务：例如，bcftools可以选择删除除一个之外的所有INFO元素。

这个工具的缺点是，它严格要求vcf格式。所以它会在你的例子中产生错误。但我认为你缩短了这篇文章的标题，你的原始文件应该没问题。为了使它对我有用，我必须调整format definition中描述的标题，在标题中为你在变体中注释的每个不同的INFO添加fileformat和INFO行：

##fileformat=VCFv4.1
##FILTER=<ID=PASS,Description="All filters passed">
##INFO=<ID=dbSNP_150,Number=0,Type=Flag,Description="has dbSNP 150 entry">
##INFO=<ID=E_Freq,Number=0,Type=Flag,Description="has E_Freq entry">
##INFO=<ID=E_1000G,Number=0,Type=Flag,Description="has E_1000G entry">
##INFO=<ID=EAS_AF,Number=A,Type=Float,Description="Allele frequency in the EAS populations">
##INFO=<ID=EUR_AF,Number=A,Type=Float,Description="Allele frequency in the EUR populations">
##INFO=<ID=AFR_AF,Number=A,Type=Float,Description="Allele frequency in the AFR populations">
##INFO=<ID=AMR_AF,Number=A,Type=Float,Description="Allele frequency in the AMR populations">
##INFO=<ID=SAS_AF,Number=A,Type=Float,Description="Allele frequency in the SAS populations">
##INFO=<ID=TSA,Number=1,Type=String,Description="No idea">
##INFO=<ID=MA,Number=1,Type=String,Description="Minor Allele">
##INFO=<ID=MAF,Number=1,Type=Float,Description="Minor Allele Frequency">
##contig=<ID=1>
##bcftools_annotateVersion=1.4-23-ga468a83+htslib-1.4-34-g8e1be4a
##bcftools_annotateCommand=annotate -x QUAL test2.vcf.gz; Date=Thu Jun 28 17:58:17 2018
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO
1       10177   rs367896724     A       AC      .       .       dbSNP_150;E_Freq;E_1000G;EAS_AF=0.3363;SAS_AF=0.4949;AFR_AF=0.4909
1       10235   rs540431307     T       TA      .       .       dbSNP_150;E_Freq;E_1000G;EAS_AF=0;AMR_AF=0.0014
1       10352   rs555500075     T       TA      .       .       dbSNP_150;EAS_AF=0.4306;EUR_AF=0.4264;SAS_AF=0.4192
1       10505   rs548419688     A       T       .       .       TSA=SNV;E_Freq;EAS_AF=0;AMR_AF=0;AFR_AF=0.0008
1       10506   rs568405545     C       G       .       .       dbSNP_150;TSA=SNV;MA=G;MAF=0.000199681;EAS_AF=0

根据您的标题，您可能需要执行两个额外步骤才能使用此工具。首先是用bgzip压缩你的输入：

bgzip phase_3.vcf

第二个是制作一个tabix-index来快速访问你的压缩文件（这会创建一个额外的文件phase_3.vcf.gz.tbi作为输出）：

tabix phase_3.vcf.gz

输入格式正确后，对bcftools的实际调用只是：

bcftools annotate -x ^INFO/EAS_AF phase_3.vcf.gz >phase_3_output.vcf

通过这些步骤，我得到的东西非常接近您想要的输出：

##fileformat=VCFv4.1
##FILTER=<ID=PASS,Description="All filters passed">
##INFO=<ID=EAS_AF,Number=A,Type=Float,Description="Allele frequency in the EAS populations">
##contig=<ID=1>
##bcftools_annotateVersion=1.4-23-ga468a83+htslib-1.4-34-g8e1be4a
##bcftools_annotateCommand=annotate -x ^INFO/EAS_AF test2.vcf.gz; Date=Thu Jun 28 17:59:04 2018
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO
1       10177   rs367896724     A       AC      .       .       EAS_AF=0.3363
1       10235   rs540431307     T       TA      .       .       EAS_AF=0
1       10352   rs555500075     T       TA      .       .       EAS_AF=0.4306
1       10505   rs548419688     A       T       .       .       EAS_AF=0
1       10506   rs568405545     C       G       .       .       EAS_AF=0

删除头部的前几行，删除带有剪切的QUAL和FILTER列，并使用sed重写0到0.0000完成您的任务：

bcftools annotate -x ^INFO/EAS_AF phase_3.vcf.gz | tail -n +7 | cut -f1-5,8 |sed 's/=0$/=0.0000/g' >phase_3_output_finished.vcf

你得到了理想的结果：

#CHROM  POS ID  REF ALT INFO
1   10177   rs367896724 A   AC  EAS_AF=0.3363
1   10235   rs540431307 T   TA  EAS_AF=0.0000
1   10352   rs555500075 T   TA  EAS_AF=0.4306
1   10505   rs548419688 A   T   EAS_AF=0.0000
1   10511   rs534229142 G   A   EAS_AF=0.0000
1   10539   rs537182016 C   A   EAS_AF=0.0000

根据您的最终目标，您甚至可以通过查找bcftools等工具选项来找到实现此步骤的方法。

因为您正在使用生物信息学数据，您可能会在biostars或bioinformatics.stackexchange等社区中找到更多帮助。

0
投票

如果不需要，摆脱第一行，然后使用pandas导入csv并命名数据框的列

    import pandas as pd
    import numpy as np
    df = pd.read_csv('/example.csv',names=['Column1','Column2'....])

编辑：这是一个python脚本，您需要在Linux系统中安装python或使用内置版本。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.