pixy软件是比较群体基因组差异的一款软件
在使用时发现只能固定窗口，不能滑动步长，这在联合比较其他算法计算的选择信号时需要用到。
所以记录一下怎么弄吧
1.思路
本身就是计算窗口内的dxy，fst这些吧，，命令可以指定窗口，下面是用法。
1
2
3
4
| usage: pixy [-h] --stats {pi,dxy,fst} [{pi,dxy,fst} ...] --vcf [VCF] --populations [POPULATIONS] [--window_size [WINDOW_SIZE]] [--bed_file [BED_FILE]] [--n_cores [N_CORES]]
[--output_folder [OUTPUT_FOLDER]] [--output_prefix [OUTPUT_PREFIX]] [--chromosomes [CHROMOSOMES]] [--interval_start [INTERVAL_START]] [--interval_end [INTERVAL_END]]
[--sites_file [SITES_FILE]] [--chunk_size [CHUNK_SIZE]] [--fst_type {wc,hudson}] [--bypass_invariant_check {yes,no}] [--version] [--citation] [--silent]
pixy: error: the following arguments are required: --stats, --vcf, --populations
|
那么手动以step空隔出来区间去计算区间内是不是也合理呢。
2.打印命令
那么多窗口手动都不现实，循环又不知道窗口合理大小，那就看一下染色体长度，打印一下命令吧。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
| #!/bin/bash
vcf_file="/home/liuxiao/1254/group1/370.vcf.gz"
window_size=150000
step_size=75000
output_intervals="intervals.txt"
chromosomes=$(bcftools view -h $vcf_file | grep "^##contig" | sed 's/##contig=<ID=//;s/,length=/\t/' | cut -f1,2)
> $output_intervals
while read -r chromosome length; do
length=$(echo $length | tr -d '>')
echo "Processing chromosome $chromosome of length $length"
start=1
end=$window_size
while [ $start -le $length ]; do
if [ $end -gt $length ]; then
end=$length
fi
echo -e "$chromosome\t$start\t$end" >> $output_intervals
start=$((start + step_size))
end=$((start + window_size - 1))
done
done <<< "$chromosomes"
echo "根据染色体规划的窗口输出: $output_intervals"
|
3.开工
知道窗口了，那这么多窗口一个一个循环计算岂不是慢死了。
parallel启动
把窗口计算完整命令打印一下
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
| #!/bin/bash
vcf_file="/home/liuxiao/1254/group1/370.vcf.gz"
populations_file="/home/liuxiao/1254/group1/vcf/dxy/pop/hu_other.txt"
output_prefix="tepedxy_outputwinstep_pophu_other"
intervals_file="intervals.txt"
n_cores=60
temp_dir=$(mktemp -d)
output_folder="$temp_dir/output"
mkdir -p $output_folder
while read -r chromosome start end; do
echo "Processing $chromosome: $start-$end"
pixy --stats dxy --vcf $vcf_file --populations $populations_file --window_size $((end-start+1)) \
--chromosomes $chromosome --interval_start $start --interval_end $end \
--n_cores $n_cores --bypass_invariant_check yes --output_folder $output_folder --output_prefix ${output_prefix}_${chromosome}_${start}_${end}
done < $intervals_file
cat $output_folder/${output_prefix}_* > ${output_prefix}_combined.txt
rm -r $temp_dir
echo "command输出到 ${output_prefix}_combined.txt"
|
如上图，根据染色体号和区间填充了我们的命令
直接bash就算了，parallel给服务器上上强度
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
| #!/bin/bash
commands_file="pixy_commands.txt"
output_prefix="tepedxy_outputwinstep_pophu_other"
output_folder="hu_other"
mkdir -p $output_folder
parallel --jobs 30 < $commands_file
wait
for chr in $(grep --only-matching --perl-regexp "(?<=--chromosomes )[^\s]+" $commands_file | sort | uniq); do
echo "Merging results for chromosome $chr"
cat ${output_prefix}_${chr}_* > ${output_folder}/combine_${chr}_combined.txt
done
wait
cat ${output_folder}/combine_*_combined.txt > ${output_folder}/${output_prefix}_final_combined.txt
wait
echo "分析完成，结果已复制到${output_folder}/${output_prefix}_final_combined.txt"
|
一小点点的不足，过程文件太多了，rm处理不了，自己分染色体rm吧
1
2
3
|
for i in {1..27}; rm *Chr$i*
rm *tmp
|
还有坑，合并的子文件标题行。。。
1
| awk 'NR==1 || !/^pop1/' combine.txt > cutcombine.txt
|
排序也不对
1
2
3
4
5
6
7
8
9
10
11
12
13
| #!/bin/bash
input_file="step.txt"
output_file="step_sorted.txt"
header_file="header.txt"
data_file="data.txt"
head -n 1 $input_file > $header_file
tail -n +2 $input_file > $data_file
sort -k3,3 -k4,4n $data_file > $output_file
cat $header_file $output_file > ${output_file}_final
echo "大功告成了，看看 ${output_file}."
|
这差不多成功了，
1
| cut -f 3,4,5,6 $output_file > clean_file
|