我有一个5000行的文本文件。每一行对应于一个唯一的FILE,全部标记为gwas*。
CHR BP FILE
chr1 12345678 gwas1
chr2 87654321 gwas2
...我有5,000个gwas*文件,它们都有唯一的文件名(如上面第4列所示- FILE),例如,gwas1如下所示:
CHR BP
chr1 12345678
chr1 12345679
chr1 12356777
...我希望从每个gwas*文件中提取所有行,其中BP值在文本文件中BP值的500,000 (上下)范围内。CHR值也必须匹配。例如:
gwas1具有CHR值chr1和BP值12,345,678。gwas1中提取所有行,这些行在CHR列中有一个chr1值,并且在11,845,678 (这个值比文本文件中的BP值下降了500,000 )和12,845,678 (从文本文件中的BP值下降了500,000 )之间有一个BP值。我可以使用以下代码手动对单个gwas文件执行此操作(但这不使用5,000行的文本文件):
export CHR="chr11"
export BP=107459522
export WINDOW=500000
awk -v CHR=$CHR -v BP_pos=$(($BP + $WINDOW)) -v BP_neg=$(($BP - $WINDOW)) 'BEGIN{FS=OFS="\t"}FNR==1 || ($1 == CHR && $2 < BP_pos && $2 > BP_neg )' gwas1 > gwas1_extract我希望能够对我的文本文件中列出的所有5000个gwas文件执行此操作。每个gwas文件的输出应该只包括(i)列CHR包含该文件的值的行(例如,对于gwas1,这是chr1),BP列包含的值在BP列文本文件中给出的值500,000以内(对于gwas1,这是12,345,678的任何一侧)。输出文件如下所示:
CHR BP
chr1 12345678
chr1 12345679
chr1 12356777awk --version = GNU 4.0.2
任何帮助都会很好!
发布于 2022-07-08 16:04:26
样本输入:
$ cat file.txt
CHR BP FILE
chr1 12345678 gwas1
chr2 87654321 gwas2
chr4 99999999 gwas4 # file gwas4 does not exist
$ cat gwas1
CHR BP
chr1 12345678 # match
chr1 12345679 # match
chr1 12356777 # match
chr1 99999999
$ cat gwas2
CHR BP
chr1 12345678
chr2 87650000 # match
$ cat gwas3 # no matches since gwas3 not in file.txt
CHR BP
chr3 2134234注意:实际文件中不存在注释
处理所有GNU awk文件的单个gwas*脚本:
awk -v diff=500000 '
function abs(x) { return (x < 0.0) ? -x : x }
BEGIN { FS=OFS="\t" }
FNR==NR { if (FNR>1) # 1st file; skip header and ...
bp_list[$3][$1]=$2 # save contents in our bp_list[FILE][CHR] array
next
}
FNR==1 { close(outfile) # close previous output file
fn=FILENAME
outfile=fn "_extract"
if (fn in bp_list) # if fn in 1st file then ...
print > outfile # print header else ...
else # skip to next input file; also addresses gwas* matching on gwas*_extract files, ie, these will be skipped, too
nextfile
next
}
fn in bp_list { if ($1 in bp_list[fn] && abs(bp_list[fn][$1] - $2) <= diff)
print > outfile
}
' file.txt gwas*注意:需要多维数组(也称为数组数组)的GNU awk
这就产生了:
$ head gwas*extract
==> gwas1_extract <==
CHR BP
chr1 12345678
chr1 12345679
chr1 12356777
==> gwas2_extract <==
CHR BP
chr2 87650000发布于 2022-07-07 12:22:42
任何帮助都会很好!
您可能会发现GNU AWK的ARGC和ARGV在您的案例中很有用,如
一个程序可以改变ARGC和ARGV的元素。每次awk到达输入文件的末尾时,它都使用ARGV的下一个元素作为下一个输入文件的名称。通过在那里存储一个不同的字符串,程序可以更改读取哪些文件。使用"-“表示标准输入。存储其他元素和递增ARGC会导致读取其他文件。
因此,您可以根据文件的内容选择要处理的文件。
https://stackoverflow.com/questions/72897318
复制相似问题