☰

「while」ループに「awk」を入れ子にして、2つのファイルを1行ずつ解析し、列の値を比較します。

tag-icon • tag-icon shell-script text-processing awk bioinformatics
「while」ループに「awk」を入れ子にして、2つのファイルを1行ずつ解析し、列の値を比較します。

awkいくつかの＆ループの組み合わせの助けが必要ですwhile。列を持つ2つの単純なファイル（通常のファイルは非常に大きい）があります。 1つはID = 10（コーディング領域（エクソン）、ここでは染色体10）の単純な間隔を表します。

#exons.bed
10  60005   60100   
10  61007   61130   
10  61200   61300   
10  61500   61650   
10  61680   61850

もう1つは、順次読み込み（=間隔はありますが小さい）、別の値を最後の列として使用することです。これは後で必要です。

#reads.bed
10  60005   60010    34 
10  61010   61020    40
10  61030   61040    22
10  61065   61070    35 
10  61100   61105    41

だから、迅速で効率的な方法で検索して、どの読み取り間隔（ファイルのどの行）とエンコード領域に属する読み取り間隔がいくつあるかを調べたいと思います。

exon 1(first interval of table 1) contains reads of line 1,2,3, etc. 
of   reads.file(2nd table)

これにより、後で各エクソンのこの行の4番目の列の値を取得できます。

各awkに対して読み取り行を1つずつ解析できないため、whileループに対するいくつかの変更が必要な場合があるコードスニペットを作成しました。ここにいる：

while read chr a b cov; do  #for the 4-column file

#if <a..b> interval of read falls inside exon interval:
awk '($2<=$a && $b <= $3) {print NR}' exons.bed >> out_lines.bed

done < reads.bed

現在a、bを手動で提供するときにawk行を実行できますが、自動的に実行したいと思います。各ペアa、bについてファイルを介して。

構文の変更や変更方法の提案をいただきありがとうございます！

フォローアップ

最後に、次のコードで問題を解決しました。

    awk 'NR==FNR{
        a[NR]=$2; 
        b[NR]=$3;
        next; }
    {  #second file
    s[i]=0; m[i]=0;  k[i]=0;              # Add sum and mean calculation
    for (i in a){                                            
       if($2>=a[i] && $3<=b[i]){         # 2,3: cols of second file here
          k[i]+=1
          print k                      #Count nb of reads found in
          out[i]=out[i]" "FNR          # keep Nb of Line of read 
          rc[i]=rc[i]" "FNR"|"$4       #keep Line and cov value of $4th col
          s[i]= s[i]+$4                #sum over coverages for each exon
          m[i]= s[i]/k[i]             #Calculate mean (k will be the No or  
                                       #reads found on i-th exon)
     }}  
    }
    END{
       for (i in out){
          print "Exon", i,": Reads with their COV:",rc[i],\
          "Sum=",s[i],"Mean=",m[i] >> "MeanCalc.txt"

    }}' exons.bed  reads.bed

出力：

   Exon 2 : Reads with their COV:  2|40 3|22 4|35 5|41 Sum= 138  Mean= 34.5
   etc.

ベストアンサー1

awk最初の問題は、そのように内部でbash変数を使用できないことです。$a内部awk評価は大地 aしかし、はaに定義されていないので空です。この問題を解決する1つの方法は、のオプションを使用して変数を定義することです。awkbashawk-v

-v var=val
--assign var=val
   Assign the value val to the variable var,  before  execution  of
   the  program  begins.  Such variable values are available to the
   BEGIN rule of an AWK program.

したがって、次のようにすることができます。

while read chr a b cov; do 
  awk -v a="$a" -v b="$b" '($2<=a && b <= $3) {print NR}' exons.bed > out$a$b 
done < reads.bed

しかし、別のエラーがあります。読み取りがエクソン内に属するには、読み取りの開始位置がエクソンの開始位置より大きく、終了位置がエクソンの終了位置より小さくなければなりません。これを使用して、$2<=a && b <= $3エクソン境界の外側から始まる読み取りを選択します。あなたが望むものです$2>=a && $3<=b。

とにかく、bashループでこれらのタスクを実行することは、各sumペアに対してa入力ファイルを一度読み取る必要があるため、非常に非効率的ですb。なぜやりませんかawk？

awk 'NR==FNR{a[NR]=$2;b[NR]=$3; next} {
        for (i in a){
           if($2>=a[i] && $3<=b[i]){
            out[i]=out[i]" "FNR 
        }}}
        END{for (i in out){
                   print "Exon",i,"contains reads of line(s)"out[i],\
                   "of reads file" 
        }}' exons.bed reads.bed

上記のスクリプトをサンプルファイルで実行すると、次の出力が生成されます。

Exon 1 contains reads of line(s) 1 of reads file
Exon 2 contains reads of line(s) 2 3 4 5 of reads file

わかりやすくするために、ここでは短縮されていない形式で同じ内容があります。

#!/usr/bin/awk -f

## While we're reading the 1st file, exons.bed
NR==FNR{
    ## Save the start position in array a and the end 
    ## in array b. The keys of the arrays are the line numbers.
    a[NR]=$2;
    b[NR]=$3; 
    ## Move to the next line, without continuing
    ## the script.
    next;
}
 ## Once we move on to the 2nd file, reads.bed
 {
     ## For each set of start and end positions
     for (i in a){
         ## If the current line's 2nd field is greater than
         ## this start position and smaller than this end position,
         ## add this line number (FNR is the current file's line number)
         ## to the list of reads for the current value of i. 
         if($2>=a[i] && $3<=b[i]){
             out[i]=out[i]" "FNR 
         }
     }
 }
 ## After both files have been processed
 END{
     ## For each exon in the out array
     for (i in out){
         ## Print the exon name and the redas it contains
         print "Exon",i,"contains reads of line(s)"out[i],
             "of reads file" 
        }

おすすめ記事