生信：种群历史有效群体大小推断 SMC++(一)
SMC++ 可以基于输入的群体 vcf 文件以及亚群样品列表进行种群历史 有效群体大小推断。
• 数据准备
1）基因组文件: genome.fasta
2）群体 vcf 文件
3）突变速率
4）群体分群信息
准备 mask 文件标记缺失数据区域
由于 smc++ 软件会将 vcf 中未出现的区域作为纯合片段处理，因此需 要准 bed 格式 mask 文件来标记数据缺失的区域。这里使用 SNPable 来处理。输入文件为基因组 fasta 文件 下面是对给定代码的详细解读，将其拆分为不同的代码块，并说明每个代码块的目标和结果：
代码块1
#!/bin/sh
source ~/.bashrc
if [ $# -lt 1 ]; then
echo "Usage: $0 <fasta> "
exit 1
fi
echo -n 'start at ';date +'%Y-%m-%d %H:%M:%S'
echo "File name is $0"
目标：检查输入参数的数量，打印开始执行的时间，并显示脚本的文件名。
结果：如果输入参数数量小于1，则输出用法提示信息并退出。打印开始执行的时间和脚本的文件名。
代码块2
mv $1 genome.fa
bwa index genome.fa
samtools faidx genome.fa
java -Xmx4G -jar /home/5k.rice/software/picard-tools-2.1.1/picard.jar CreateSequenceDictionary R=genome.fa O=genome.dict
目标：将输入的fasta文件重命名为
genome.fa，使用bwa对genome.fa创建索引，使用samtools对genome.fa创建索引，使用picard工具创建genome.fa的序列字典。结果：将输入的fasta文件重命名为
genome.fa，创建了genome.fa的索引文件和序列字典文件。
代码块3
splitfa genome.fa 35 | split -l 20000000 -d - read.split
目标：将
genome.fa按照35个字符为一行进行分割，并按照每个文件包含的行数（20000000行）进行分割。结果：生成多个以
read.split开头的文件，每个文件包含20000000行以及35个字符为一行的片段。
代码块4
ls ./read.split*[0-9] | while read aa
do
echo "#!/bin/sh
source ~/.bashrc
bwa aln -t 4 -R 1000000 -O 3 -E 3 genome.fa $aa > $aa.sai
bwa samse genome.fa $aa.sai $aa |gzip > $aa.sam.gz
wait && echo "koflagsamgz is done"
" > $aa.bwa.sh
done
目标：为每个
read.split开头的文件生成对应的bwa.sh脚本。结果：生成多个以
read.split开头的bwa.sh脚本文件，每个脚本文件包含bwa aln和bwa samse命令的执行，并将结果压缩保存到.sam.gz文件中。
代码块5
ls read.split*.bwa.sh | while read file
do
sbatch $file 1>$file.log 2>$file.err
done
echo "sbatch is done"
目标：遍历所有的
bwa.sh脚本文件，并通过sbatch命令提交任务。结果：使用
sbatch命令提交每个bwa.sh脚本文件作为一个任务，输出任务日志和错误文件。
代码块6
num=`ls read.split*.bwa.sh|wc -l`
while true
do
count=$(grep "koflagsamgz" *out | wc -l)
if [ $count -eq $num ]
then
break
else
sleep 60
fi
done
echo "Continue executing the commands."
目标：等待所有任务完成。
结果：持续检查输出文件中是否存在字符串"koflagsamgz"，直到其出现的次数与任务数量相等时，跳出循环。然后输出"Continue executing the commands."。
代码块7
gzip -dc read.split*.sam.gz | gen_raw_mask.pl > rawMask_35.fa
## 设置过滤阈值r=0.5
gen_mask -l 35 -r 0.5 rawMask_35.fa > mask_35_50.fa
目标：从压缩的
.sam.gz文件中解压缩并
生成rawMask_35.fa文件，然后使用gen_mask工具生成mask_35_50.fa文件。
结果：生成了
rawMask_35.fa和mask_35_50.fa两个文件。
代码块8
apply_mask_s mask_35_50.fa genome.fa > genome.mask.fa
目标：将
mask_35_50.fa中的掩码信息应用到genome.fa，生成带有掩码的genome.mask.fa文件。结果：生成了带有掩码信息的
genome.mask.fa文件。
代码块9
cp /home/perl5/perl/get_lcase.pl ./
perl get_lcase.pl genome.mask.fa > genome.mask.bed
目标：复制
/home/perl5/perl/get_lcase.pl文件到当前目录，并使用perl脚本处理genome.mask.fa文件，生成genome.mask.bed文件。结果：复制了
get_lcase.pl文件到当前目录，并通过该脚本处理genome.mask.fa文件，生成了genome.mask.bed文件。
代码块10
bgzip genome.mask.bed
tabix genome.mask.bed.gz
目标：将
genome.mask.bed文件进行压缩，并生成索引文件。结果：生成了压缩的
genome.mask.bed.gz文件，并生成了对应的索引文件。
代码块11
wait && echo "koflags is done"
## clear adm
rm -rf read.split*
rm -rf *err
rm -rf *out
echo -n 'end at ';date +'%Y-%m-%d %H:%M:%S'
目标：等待所有任务完成，清理文件，并打印结束执行的时间。
结果：等待所有任务完成，删除
read.split*文件和错误、输出文件，最后打印结束执行的时间。
国靓的文案
主要分享：1：基于R语言、python和shell的数据分析，数据可视化及生物信息分析；2：零基础学英语从How are you开始学；3：运动健康。（能量是守恒的，喜欢是互相的，关注我，世界上就多了一个爱你的人！）
