使用 linux 或 python 从文件中提取特定的列和字符串
Extract the specific columns and strings from a file using linux or python
我在处理 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
如果不需要,请删除第一行,然后使用 pandas 导入 csv 并命名数据框的列
import pandas as pd
import numpy as np
df = pd.read_csv('/example.csv',names=['Column1','Column2'....])
编辑:这是一个 python 脚本,您需要在 linux 系统中安装 python 或使用内置版本。
您的文件在一个非常普通的 bioinformatics-format(vcf) 中。因此,生物信息学中有专门为解决您的任务而设计的工具:
例如 bcftools 有一个选项可以删除除一个之外的所有 INFO-elements。
此工具的缺点是,它严格要求 vcf-format。所以它会在你的例子中产生错误。但我认为你为这个 post 缩短了 header 并且它在你的原始文件上应该没问题。为了使它对我可用,我必须像 format definition 中描述的那样调整 header,方法是在 header 中为您在变体中注释的每个不同信息添加文件格式和 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
根据您的 header,您可能需要执行两个额外的步骤才能使用此工具。首先是使用 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
使用 head 删除前几行,使用 cut 删除 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 or bioinformatics.stackexchange
等相关社区中找到更多帮助
我在处理 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
如果不需要,请删除第一行,然后使用 pandas 导入 csv 并命名数据框的列
import pandas as pd
import numpy as np
df = pd.read_csv('/example.csv',names=['Column1','Column2'....])
编辑:这是一个 python 脚本,您需要在 linux 系统中安装 python 或使用内置版本。
您的文件在一个非常普通的 bioinformatics-format(vcf) 中。因此,生物信息学中有专门为解决您的任务而设计的工具: 例如 bcftools 有一个选项可以删除除一个之外的所有 INFO-elements。
此工具的缺点是,它严格要求 vcf-format。所以它会在你的例子中产生错误。但我认为你为这个 post 缩短了 header 并且它在你的原始文件上应该没问题。为了使它对我可用,我必须像 format definition 中描述的那样调整 header,方法是在 header 中为您在变体中注释的每个不同信息添加文件格式和 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
根据您的 header,您可能需要执行两个额外的步骤才能使用此工具。首先是使用 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
使用 head 删除前几行,使用 cut 删除 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 or bioinformatics.stackexchange
等相关社区中找到更多帮助