我想将8个文件(每个文件代表一个染色体)分解成大约5个4e8行的块,每个文件大约有2e9行。这些是VCF文件(https://en.wikipedia.org/wiki/Variant_Call_Format),它有一个标题,然后是遗传变异,所以我需要保留每个文件的标题,并将它们重新附加到染色体特定的标题。我是在HPC上用linux来做这件事的。
在使用之前,我已经用一个文件完成了这项工作:
#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循环可以更好地完成这项工作,但是,语法对我来说有点混乱。
我试过了:
#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这会将所有内容放入相同的标题中。如何将输出放入每个染色体的唯一标题中?类似地,这将如何在分裂变体并将头部重新连接到每条染色体的每个块上发挥作用?
期望的输出将是:
chr001_chunk1.vcf
chr001_chunk2.vcf
chr001_chunk3.vcf
chr001_chunk4.vcf
chr001_chunk5.vcf
...
chr008_chunk5.vcf 其中每个vcf块具有来自它们各自的染色体“双亲”的头部。
非常感谢
发布于 2020-08-18 22:28:32
#!/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 0https://stackoverflow.com/questions/63465595
复制相似问题