【问题标题】:Bash Awk Select lines from a file with columns in common present in several other filesBash Awk 从一个文件中选择行,这些行在其他几个文件中存在共同的列
【发布时间】:2021-02-02 04:52:17
【问题描述】:

我有一个包含制表符分隔符 .tsv 文件的文件夹,如下所示:

Sample_1.tsv 文件:

Gene    Center  Start   End Strand  Ref Sample  ID_Sample   HGVps   Other_Columns
AAA .   111111  111111  +   T   C   Test    p.A123B NA
AAA .   111112  111112  +   C   A   Test    p.C456D NA
BBB .   222222  222222  +   A   T   Test    p.E789F NA
CCC .   333331  333331  +   G   C   Test    p.G10H  NA
CCC .   333332  333332  +   A   T   Test    p.I11J  NA
CCC .   333333  333333  +   T   C   Test    p.K12L  NA

我在另一个文件夹中有每个基因的几个文件(后来称为 Genes.tsv)。

AAA 基因的 AAA.tsv 文件:

Coordinates Some_col_1  Change  Some_col_2  Consequence Other_Columns
chr1:111111-111111  NA  chr1:g.111111T>C    NA  Ms AAA A123B    NA

BBB 基因的 BBB.tsv 文件:

Coordinates Some_col_1  Change  Some_col_2  Consequence Other_Columns
chr2:222222-222222  NA  chr2:g.222222A>T    NA  Syn BBB E789F   NA

CCC 基因的 CCC.tsv 文件:

Coordinates Some_col_1  Change  Some_col_2  Consequence Other_Columns
chr3:333332-333332  NA  chr3:g.333332A>T    NA  Nns CCC I11J    NA
chr3:333339-333339  NA  chr3:g.333339T>C    NA  Syn CCC K12L    NA

等等。对于其他基因。

请注意,在“后果”列中,三个元素由空格分隔,包括基因名称和我感兴趣的 HGVps 代码。

我将这些简化的表格作为示例放在这里,但实际上我有更多的列。请参阅下图了解正确的列号:

Image-1:Sample_1.tsv_with_correct_column_numbers注意:黄色是我想要保留的行。)

Image-2:Genes.tsv_with_correct_column_numbers注意:绿色是与 Sample_1 相同的元素。)

基本上,我想将具有其他 Genes.tsv 文件中的列的 Sample_1.tsv 的所有行保存在一个新文件中。标准是 Genes Columns+Start+End+Ref+Sample 相同或 仅 HVps 列(如果其他列不匹配)。

我想要获得的东西:

Gene    Center  Start   End Strand  Ref Sample  ID_Sample   HGVps   Other_Columns
AAA .   111111  111111  +   T   C   Test    p.A123B NA
BBB .   222222  222222  +   A   T   Test    p.E789F NA
CCC .   333332  333332  +   A   T   Test    p.I11J  NA
CCC .   333333  333333  +   T   C   Test    p.K12L  NA

我认为可以使用 bash/awk/grep 来选择感兴趣的行。

我不完全知道如何做到这一点,但我想继续的方式应该是这样的:

Create a new file called Sample_1_ok.tsv
Add column names in Sample_1_ok.tsv

If Gene+Start+End+Ref+Sample columns from Sample_1.tsv are the same in Genes.tsv files:
Append Sample_1.tsv line in Sample_1_ok.tsv

Else if HGVps column from Sample_1.tsv is the same in Genes.tsv files:
Append Sample_1.tsv line in Sample_1_ok.tsv

您对如何进行有什么建议吗?

非常感谢!

【问题讨论】:

  • 您要处理多少个<gene>.tsv 文件?您能否提供一些有关哪些值与 <gene>.tsv 文件中的哪些列名匹配的详细信息(例如,AAA.tsv 中有 3 个 111111 副本 ...这是 StartEnd) ?您是否测试过任何代码,如果可以,请分享一下?第一个想法是awk ...将<gene>.fsv viles 解析为一个数组,然后处理Sample_1.fsv 以在数组中匹配;如果“太多”<gene>.fsv 文件,则可能将<gene>.fsv 数据重新格式化为与Sample_1.fsv 格式相同的单个文件,然后comm 数据
  • 我有 ~ 50 个<gene>.tsv 文件。对于<gene>.tsv 文件的列的详细信息,Coordinates 列对应于chr#:Start-End coordinates。因为它是单字母突变,StartEnd 通常是相同的(但并非总是如此,这就是我想区分这两个坐标的原因)。 Change 列,我只想获取 > 符号前后的两个大写字母。这两个字母中的第一个对应于Ref,第二个对应于Sample。至少Consequence 列包含Gene 的名称(AAA、BBB、CCC、...)和HGVps
  • @markp-fuso 通过使用 awk 匹配,我从 <gene>.tsv 中提取了感兴趣的列:awk -F'\t' -v OFS='\t' 'match($1, /([0-9]+)\-([0-9]+)/, a) && match($3, /(.)>(.)/, b) && match($5, /([A-Z]+) ([A-Z][0-9]+[A-Z])/, c) {print c[1],a[1],a[2],b[1],b[2],c[2]}' <gene>.tsv > <gene_filtered>.tsv 所以下一步是比较 Sample_1.tsv 和这些新表?
  • fwiw,我建议将您的 awk 代码移到问题中;让其他人更容易理解这个问题,如果他们不必去挖掘 cmets 试图将所有东西拼凑在一起
  • @markp-fuso 好的,我会记住您对代码的推荐。

标签: bash csv awk compare


【解决方案1】:

一个(有点冗长)awk 解决方案:

awk '
/Consequence/   { gfmt=1 ; next }                          # header record of a "<gene>.tsv" file so flag as needing to parse based on a g(ene) file format; skip to next line of input
/ID_Sample/     { gfmt=0 }                                 # header record of "Sample_1.tsv" file so clear g(ene) file format flag; continue to process this line of input

# gfmt==1 => parse data based on the g(ene) file format

gfmt            { gene=$6                                  # parse out gene name

                  split($1,a,/[:-]/)                       # parse out start and end values
                  gstart=a[2]
                  gend=a[3]

                  split($3,a,">")                          # parse out ref and sample values
                  gref=substr(a[1],length(a[1]),1)         # assumes ref is always a single character
                  gsample=a[2]

                  hgvps="p."$7                             # parse out HGVps and prefix with "p."; assumes prefix is always "p."

                  garray[gene,gstart,gend,gref,gsample]    # use fields as multidimensional index; to be used for testing a match based on gene data
                  hgarray[gene,hgvps]                      # use fields as multidimentional index; to be used for testing a match based on HGVps
                }

# gfmt==0 => parse data based on Sample_1.tsv file format

! gfmt          { if ( FNR==1 ) { print $0 ; next }        # print header record; skip to next line of input

                  gene=$1                                  # parse out gene values;
                  gstart=$3                                # could use field numbers in the
                  gend=$4                                  # follow-on "if" logic but wanted
                  gref=$6                                  # to give names to fields
                  gsample=$7                               # as form of documentation 

                  hgvps=$9

                  # if we have an index match in either array then print the current line to stdout

                  if ( (gene,gstart,gend,gref,gsample) in garray || (gene,hgvps) in hgarray )
                      { print $0 } 
                }

' *.tsv xx/Sample_1.tsv

注意事项

  • 移除 cmets 以整理代码。
  • OP 可以换入 awk/match 代码以解析出 &lt;gene&gt;.tsv
  • 上面的代码是从&lt;gene&gt;.tsv文件所在的目录运行的
  • 上面的代码假设所有&lt;gene&gt;.tsv文件都被处理(因此'*.tsv'
  • Sample_1.tsv 文件是最后一个提供给awk 的文件;在这种情况下,我的 Sample_1.tsv 文件位于单独的子目录中

AAA.tsvBBB.tsvCCC.tsvxx/Sample_1.tsv 运行上述代码会生成以下输出:

Gene    Center  Start   End Strand  Ref Sample  ID_Sample   HGVps   Other_Columns
AAA .   111111  111111  +   T   C   Test    p.A123B NA
BBB .   222222  222222  +   A   T   Test    p.E789F NA
CCC .   333332  333332  +   A   T   Test    p.I11J  NA
CCC .   333333  333333  +   T   C   Test    p.K12L  NA

【讨论】:

    【解决方案2】:

    感谢您的回答并向我展示了一种比较具有不同结构的两个文件的有趣方法。它非常适合我的示例表。

    正如您所提到的,我已经更改了awk/match 代码中的列名,以便区分&lt;gene&gt;.tsvSample_1.tsv 文件。它奏效了。

    我在真实数据中遇到的代码中困难的部分是选择基因名称和 hgvps 代码,它们存在于同一个单元格中并用空格分隔。我更喜欢使用补充 awk match 表达式来选择它们。我还在开头添加了awk -F'\t'paramater。

    这是我对脚本的修改(全部都在第一部分中):

    awk -F'\t' '
    /Consequence/   { gfmt=1 ; next }
    /ID_Sample/     { gfmt=0 }
    
    # gfmt==1 => parse data based on the g(ene) file format
    
    gfmt==1           { match($5, / ([A-Z0-9]{2,})/, a)
                      gene=a[1]
                      print gene
    
                      split($1,a,/[:-]/)
                      gstart=a[2]
                      gend=a[3]
    
                      split($3,a,">")
                      gref=substr(a[1],length(a[1]),1)
                      gsample=a[2]
    
                      match($5, /([A-Z][0-9]+[A-Z])/, a)
                      hgvps="p," a[1]
    
                      garray[gene,gstart,gend,gref,gsample]
                      hgarray[gene,hgvps]
    
                    }
    # gfmt==0 => parse data based on Sample_1.tsv file format
    
    gfmt==0           { if ( FNR==1 ) { print $0 ; next }
    
                      gene=$1
                      gstart=$3
                      gend=$4
                      gref=$6
                      gsample=$7
    
                      hgvps=$9
                      if ( (gene,gstart,gend,gref,gsample) in garray || (gene,hgvps) in hgarray )
                          { print $0 }
                    }
    
    ' *.tsv xx/Sample_1.tsv
    
    

    再次感谢您的帮助!

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2011-01-16
      • 1970-01-01
      • 1970-01-01
      • 2023-03-10
      • 2011-10-18
      • 1970-01-01
      相关资源
      最近更新 更多