首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >在具有不同输出文件的多个文件上运行数组linux

在具有不同输出文件的多个文件上运行数组linux
EN

Stack Overflow用户
提问于 2020-08-18 17:20:55
回答 1查看 48关注 0票数 0

我想将8个文件(每个文件代表一个染色体)分解成大约5个4e8行的块,每个文件大约有2e9行。这些是VCF文件(https://en.wikipedia.org/wiki/Variant_Call_Format),它有一个标题,然后是遗传变异,所以我需要保留每个文件的标题,并将它们重新附加到染色体特定的标题。我是在HPC上用linux来做这件事的。

在使用之前,我已经用一个文件完成了这项工作:

代码语言:javascript
复制
#grab the header
head -n 10000 my.vcf | grep "^#" >header
#grab the non header lines
grep -v "^#" my.vcf >variants
#split into chunks with 40000000 lines
split -l 40000000 variants
#reattach the header to each and clean up
for i in x*;do cat header $i >$i.vcf && rm -f $i;done
rm -f header variants

我可以手动处理所有8条染色体,但我在具有数组功能的HPC中工作,我觉得使用for循环可以更好地完成这项工作,但是,语法对我来说有点混乱。

我试过了:

代码语言:javascript
复制
#filelist is a list of the 8 chromosome files i.e. chr001.vcf, chr002.vcf...chr0008.vcf 
for f in 'cat filelist.txt'; head -n 10000 my.vcf | grep "^#" >header; done

这会将所有内容放入相同的标题中。如何将输出放入每个染色体的唯一标题中?类似地,这将如何在分裂变体并将头部重新连接到每条染色体的每个块上发挥作用?

期望的输出将是:

代码语言:javascript
复制
chr001_chunk1.vcf
chr001_chunk2.vcf
chr001_chunk3.vcf
chr001_chunk4.vcf
chr001_chunk5.vcf
...
chr008_chunk5.vcf 

其中每个vcf块具有来自它们各自的染色体“双亲”的头部。

非常感谢

EN

回答 1

Stack Overflow用户

回答已采纳

发布于 2020-08-18 22:28:32

代码语言:javascript
复制
#!/bin/bash

#
# scan the current directory for chr[0-9]*.vcf
# extract header lines (^#)
# extract variants (non-header lines) and split to 40m partial files
# combine header with each partial file
#

# for tuning
lines=40000000

vcf_list=(chr[0-9]*.vcf)
if [ ${#vcf_list} -eq 0 ]; then
    echo no .vcf files
    exit 1
fi

tmpv=variants
hdr=header

for chrfile in "${vcf_list[@]}"; do
    # isolate without . extn
    base=${chrfile%%.*}
    echo $chrfile

    # extract header lines
    head -1000 $chrfile | grep "^#" > $hdr

    # extract variants
    grep -v "^#" $chrfile > $tmpv

    #
    # split variants into files with max $lines;
    # output files are created with a filter to combine header data and
    # partial variant data in 1 pass, avoiding additional file I/O;
    # output files are named with a leading 'p' to support multiple
    # runs without filename collision
    #
    split -d -l $lines $tmpv p${base}_chunk --additional-suffix=.vcf \
        --filter="cat $hdr - > \$FILE; echo \"  \$FILE\""
done

rm -f $tmpv $hdr

exit 0
票数 2
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/63465595

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档