主页 / unix / 问题 / 427831
Accepted
Age87
Asked: 2018-03-03 18:38:11 +0800 CST

使用文件子集中的信息从文件中提取单词（多个步骤）

  • 772

我有两个大文件，如下所示：

文件 1：

NW_006502347.1 316684
NW_006527876.1 351
NW_006502151.1 27628
NW_006526579.1 232
NW_006525259.1 132
NW_006501641.1 437014
NW_006525259.1 378
NW_006523082.1 215
NW_006522424.1 153
NW_006522101.1 815
NW_006521985.1 505
NW_006521985.1 527
NW_006521722.1 920
NW_006521525.1 73
NW_006521432.1 258
NW_006521302.1 938
NW_006521272.1 585
NW_006521272.1 745
NW_006521038.1 202
NW_006519846.1 1528
NW_006519837.1 10215

文件 2：

NW_006502347.1  Gnomon  CDS   305319  305340   .  +  0  ID=cds43608  Parent=rna48098  Dbxref=GeneID:102908761,Genbank:XP_006997436.2  Name=XP_006997436.2  gbkey=CDS  gene=LOC102908761  partial=true  product=histone deacetylase 4-like  protein_id=XP_006997436.2
NW_006501037.1  Gnomon  gene  6936    115174   .  -  .  ID=gene0                      Dbxref=GeneID:102922816  Name=Efl1  gbkey=Gene  gene=Efl1  gene_biotype=protein_coding  partial=true  start_range=.,6936
NW_006501037.1  Gnomon  mRNA  6936    115174   .  -  .  ID=rna0      Parent=gene0     Dbxref=GeneID:102922816,Genbank:XM_006970114.2  Name=XM_006970114.2  gbkey=mRNA  gene=Efl1  model_evidence=Supporting evidence includes similarity to: 5 mRNAs%2C 7 Proteins%2C and 99%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 16 samples with support for all annotated introns  partial=true  product=elongation factor like GTPase 1  start_range=.,6936  transcript_id=XM_006970114.2
NW_006501037.1  Gnomon  exon  115095  115174   .  -  .  ID=id1       Parent=rna0      Dbxref=GeneID:102922816,Genbank:XM_006970114.2  gbkey=mRNA  gene=Efl1  partial=true  product=elongation factor like GTPase 1  transcript_id=XM_006970114.2
NW_006501037.1  Gnomon  exon  114246  114355   .  -  .  ID=id2       Parent=rna0      Dbxref=GeneID:102922816,Genbank:XM_006970114.2  gbkey=mRNA  gene=Efl1  partial=true  product=elongation factor like GTPase 1  transcript_id=XM_006970114.2
NW_006502347.1  Gnomon  mRNA  272091  477077   .  +  .  ID=rna48098  Parent=gene26399 Dbxref=GeneID:102908761,Genbank:XM_006997374.2  Name=XM_006997374.2  end_range=477077,.  gbkey=mRNA  gene=LOC102908761  model_evidence=Supporting evidence includes similarity to: 1 mRNA%2C and 90%25 coverage of the annotated genomic feature by RNAseq alignments  partial=true  product=histone deacetylase 4-like  transcript_id=XM_006997374.2
NW_006501037.1  Gnomon  exon  92339   92472    .  -  .  ID=id5       Parent=rna0      Dbxref=GeneID:102922816,Genbank:XM_006970114.2  gbkey=mRNA  gene=Efl1  partial=true  product=elongation factor like GTPase 1  transcript_id=XM_006970114.2
NW_006501037.1  Gnomon  exon  90969   91106    .  -  .  ID=id6       Parent=rna0      Dbxref=GeneID:102922816,Genbank:XM_006970114.2  gbkey=mRNA  gene=Efl1  partial=true  product=elongation factor like GTPase 1  transcript_id=XM_006970114.2
NW_006501037.1  Gnomon  exon  89261   89475    .  -  .  ID=id7       Parent=rna0      Dbxref=GeneID:102922816,Genbank:XM_006970114.2  gbkey=mRNA  gene=Efl1  partial=true  product=elongation factor like GTPase 1  transcript_id=XM_006970114.2
NW_006502151.1  Gnomon  exon  26099   27657    .  -  .  ID=id586652  Parent=rna47002  Dbxref=GeneID:102918658,Genbank:XM_006996663.2  gbkey=mRNA  gene=Rftn1  product=raftlin%2C lipid raft linker 1  transcript_id=XM_006996663.2
NW_006501641.1  Gnomon  mRNA  393496  438556   .  +  .  ID=rna40001  Parent=gene21212 Dbxref=GeneID:102913870,Genbank:XM_015986269.1  Name=XM_015986269.1  gbkey=mRNA  gene=LOC102913870  model_evidence=Supporting evidence includes similarity to: 9 mRNAs%2C 5 Proteins%2C and 81%25 coverage of the annotated genomic feature by RNAseq alignments  product=transmembrane protein 189  transcript_id=XM_015986269.1
NW_006501053.1  Gnomon  exon  5104713 5104872  .  +  .  ID=id45206   Parent=rna3590   Dbxref=GeneID:102916769,Genbank:XR_001580019.1  gbkey=misc_RNA  gene=Mycbpap  product=MYCBP associated protein%2C transcript variant X4  transcript_id=XR_001580019.1
NW_006501053.1  Gnomon  exon  5104959 5105062  .  +  .  ID=id45207   Parent=rna3590   Dbxref=GeneID:102916769,Genbank:XR_001580019.1  gbkey=misc_RNA  gene=Mycbpap  product=MYCBP associated protein%2C transcript variant X4  transcript_id=XR_001580019.1
NW_006501053.1  Gnomon  exon  5105698 5105881  .  +  .  ID=id45208   Parent=rna3590   Dbxref=GeneID:102916769,Genbank:XR_001580019.1  gbkey=misc_RNA  gene=Mycbpap  product=MYCBP associated protein%2C transcript variant X4  transcript_id=XR_001580019.1
NW_006501053.1  Gnomon  exon  5106131 5106246  .  +  .  ID=id45209   Parent=rna3590   Dbxref=GeneID:102916769,Genbank:XR_001580019.1  gbkey=misc_RNA  gene=Mycbpap  product=MYCBP associated protein%2C transcript variant X4  transcript_id=XR_001580019.1

我想使用文件 1 中的信息来提取第 13 或 14 列下“gene=”后面的单词（例如“Efl1”）。更具体地说，我想：

步骤 1)比较文件 1 中第 1 列的标签（例如 NW_006527876.1）与文件 2 中第 1 列的标签，提取第 1 列（文件 2）与文件 1 的第 1 列匹配的所有行。

如您所见，第 1 列（文件 2）的标签重复，因此文件 1 的每个标签都会有多个匹配项。

文件 2 的第 4 列和第 5 列代表一个区间，第 4 列是区间的开始，第 5 列是区间的结束。文件 1 的第 2 列表示这些间隔之间的数字。

步骤 2) 从步骤 1 中分离出来的行中，提取第 2 列（文件 1）下的编号位于文件 2 的第 4 列和第 5 列表示的区间之间的行。

这远远超出了我所知道的范围，但下面是关于命令可能是什么样子的想法。

awk '{ print $1 }' file 1 |  
awk `$4 *(file2)* < $2 *(file1)*' | awk '$5 *(file2)* > $2 *(file1)*' > output.tsv

输出应该在第 1 列下具有具有唯一标签的行。

第 3 步）从上面创建的 output.tsv 文件中，我想提取第 13 列或第 14 列（见下文）下“gene=”中等号后面的单词，这样我就得到了一个文件用等号后面的词。

最终输出文件（基于此示例）：

LOC102908761
Rftn1
LOC102913870

回答：

while read -r id pos; do awk -v id="$id" -v pos="$pos" '$1 == id && pos > $4 && pos < $5 { if (gensub(/.*gene=([A-Za-z0-9]*).*/, "\\1", 1) !~ /\s/) print gensub(/.*gene=([A-Za-z0-9]*).*/, "\\1", 1); }' < file2.txt; done <file1.txt > gene_hits.txt

接下来要去除重复（多个间隔中的相同基因座）同时保持不同基因座映射到同一基因，请执行以下操作：

perl -ne 'print if ++$k{$_}==1' A_gene_hits.txt > A_genes.txt
awk sed

1 个回答

  1. Best Answer

    假设您有空格作为分隔符：

    $ while read -r id pos; do awk -v id="$id" -v pos="$pos" '$1 == id && pos > $4 && pos < $5 { print gensub(/.*gene=([A-Za-z0-9]*).*/, "\\1", 1) }' <file2; done <file1
LOC102908761
Rftn1
LOC102913870

    解释

    • while read -r id pos; do FOO; done <file1file1：这将逐行读取，并将第一个字段（例如NW_006502347.1）放入 shell 变量$id中，将第二个字段（例如316684）放入 shell 变量$pos中。然后它FOO为每一行运行。
    • awk -v id="$id" -v pos="$pos" 'BAR' <file2：对于每一行file1，我们将运行一个将运行的awk命令BAR。这将搜索file2匹配的部分。我们需要告诉这个awk脚本来自 shell 的两个“外部”变量。即 awk 变量id被赋予与 shell 变量相同的值$id，而 awk 变量pos和 shell 变量也被赋予相同的值$pos。
    • $1 == id && pos > $4 && pos < $5：这是awk脚本的“条件”部分。如果满足这些条件，则将运行以下命令。在这里，我们检查 的第一个字段是否与当前行的$1相同file2，并且是否在 的第 4和第 5 个字段之间。idfile1pos$4$5file2
    • { print gensub(/.*gene=([A-Za-z0-9]*).*/, "\\1", 1) }：如果满足上述条件，则此代码将运行。我们想用gensubfirst 进行替换。这将搜索gene=后跟任意长度的字母数字字符串([A-Za-z0-9]*)。此字母数字字符串由括号(捕获。我们还将“搜索”完整字符串)前后的所有字符。因此，这会“搜索”整行，将其替换为（第一个也是唯一的）捕获组，即 之后的字母数字字符串。替换第一次出现的最终方法，尽管这没有多大意义，因为我假设每行只有一个匹配项。.*gene=([A-Za-z0-9]*)"\\1"gene=1gene=

    制表符分隔的版本

    一般来说，我更喜欢使用制表符分隔的文件，尤其是对于我认为是 GFF/GTF 文件的文件。这允许区分空间，尤其是在第 9 个字段中。

    while IFS=$'\t' read -r id pos; do awk -F'\t' -v id="$id" -v pos="$pos" '$1 == id && pos > $4 && pos < $5 { print gensub(/.*gene=([A-Za-z0-9]*).*/, "\\1", 1) }' <file2.tsv ; done <file1.tsv

    对脚本的修改是显式拆分选项卡上的 shell 行IFS=$'\t'，并awk用-F'\t'.

    • 2

相关问题