SMC++ 可以基于輸入的群體 vcf 文件以及亞群樣品列表進(jìn)行種群歷史 有效群體大小推斷。?

? 數(shù)據(jù)準(zhǔn)備?

1）基因組文件: genome.fasta?

2）群體 vcf 文件

3）突變速率?

4）群體分群信息

準(zhǔn)備 mask 文件標(biāo)記缺失數(shù)據(jù)區(qū)域

由于 smc++ 軟件會(huì)將 vcf 中未出現(xiàn)的區(qū)域作為純合片段處理，因此需 要準(zhǔn) bed 格式 mask 文件來標(biāo)記數(shù)據(jù)缺失的區(qū)域。這里使用 SNPable 來處理。輸入文件為基因組 fasta 文件 下面是對(duì)給定代碼的詳細(xì)解讀，將其拆分為不同的代碼塊，并說明每個(gè)代碼塊的目標(biāo)和結(jié)果：

代碼塊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"


• 目標(biāo)：檢查輸入?yún)?shù)的數(shù)量，打印開始執(zhí)行的時(shí)間，并顯示腳本的文件名。

• 結(jié)果：如果輸入?yún)?shù)數(shù)量小于1，則輸出用法提示信息并退出。打印開始執(zhí)行的時(shí)間和腳本的文件名。

代碼塊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


• 目標(biāo)：將輸入的fasta文件重命名為genome.fa，使用bwa對(duì)genome.fa創(chuàng)建索引，使用samtools對(duì)genome.fa創(chuàng)建索引，使用picard工具創(chuàng)建genome.fa的序列字典。

• 結(jié)果：將輸入的fasta文件重命名為genome.fa，創(chuàng)建了genome.fa的索引文件和序列字典文件。

代碼塊3

splitfa?genome.fa?35?|?split?-l?20000000?-d?-?read.split


• 目標(biāo)：將genome.fa按照35個(gè)字符為一行進(jìn)行分割，并按照每個(gè)文件包含的行數(shù)（20000000行）進(jìn)行分割。

• 結(jié)果：生成多個(gè)以read.split開頭的文件，每個(gè)文件包含20000000行以及35個(gè)字符為一行的片段。

代碼塊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


• 目標(biāo)：為每個(gè)read.split開頭的文件生成對(duì)應(yīng)的bwa.sh腳本。

• 結(jié)果：生成多個(gè)以read.split開頭的bwa.sh腳本文件，每個(gè)腳本文件包含bwa aln和bwa samse命令的執(zhí)行，并將結(jié)果壓縮保存到.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"


• 目標(biāo)：遍歷所有的bwa.sh腳本文件，并通過sbatch命令提交任務(wù)。

• 結(jié)果：使用sbatch命令提交每個(gè)bwa.sh腳本文件作為一個(gè)任務(wù)，輸出任務(wù)日志和錯(cuò)誤文件。

代碼塊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."


• 目標(biāo)：等待所有任務(wù)完成。

• 結(jié)果：持續(xù)檢查輸出文件中是否存在字符串"koflagsamgz"，直到其出現(xiàn)的次數(shù)與任務(wù)數(shù)量相等時(shí)，跳出循環(huán)。然后輸出"Continue executing the commands."。

代碼塊7

gzip?-dc?read.split*.sam.gz?|?gen_raw_mask.pl?>?rawMask_35.fa
##?設(shè)置過濾閾值r=0.5
gen_mask?-l?35?-r?0.5?rawMask_35.fa?>?mask_35_50.fa


• 目標(biāo)：從壓縮的.sam.gz文件中解壓縮并

生成rawMask_35.fa文件，然后使用gen_mask工具生成mask_35_50.fa文件。

• 結(jié)果：生成了rawMask_35.fa和mask_35_50.fa兩個(gè)文件。

代碼塊8

apply_mask_s?mask_35_50.fa?genome.fa?>?genome.mask.fa


• 目標(biāo)：將mask_35_50.fa中的掩碼信息應(yīng)用到genome.fa，生成帶有掩碼的genome.mask.fa文件。

• 結(jié)果：生成了帶有掩碼信息的genome.mask.fa文件。

代碼塊9

cp?/home/perl5/perl/get_lcase.pl?./
perl?get_lcase.pl?genome.mask.fa?>?genome.mask.bed


• 目標(biāo)：復(fù)制/home/perl5/perl/get_lcase.pl文件到當(dāng)前目錄，并使用perl腳本處理genome.mask.fa文件，生成genome.mask.bed文件。

• 結(jié)果：復(fù)制了get_lcase.pl文件到當(dāng)前目錄，并通過該腳本處理genome.mask.fa文件，生成了genome.mask.bed文件。

代碼塊10

bgzip?genome.mask.bed
tabix?genome.mask.bed.gz


• 目標(biāo)：將genome.mask.bed文件進(jìn)行壓縮，并生成索引文件。

• 結(jié)果：生成了壓縮的genome.mask.bed.gz文件，并生成了對(duì)應(yīng)的索引文件。

代碼塊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'


• 目標(biāo)：等待所有任務(wù)完成，清理文件，并打印結(jié)束執(zhí)行的時(shí)間。

• 結(jié)果：等待所有任務(wù)完成，刪除read.split*文件和錯(cuò)誤、輸出文件，最后打印結(jié)束執(zhí)行的時(shí)間。
  
  

國靚的文案

主要分享：1：基于R語言、python和shell的數(shù)據(jù)分析，數(shù)據(jù)可視化及生物信息分析；2：零基礎(chǔ)學(xué)英語從How are you開始學(xué)；3：運(yùn)動(dòng)健康。（能量是守恒的，喜歡是互相的，關(guān)注我，世界上就多了一個(gè)愛你的人?。?/span>