03岳颇、測序數(shù)據(jù)批量比對到參考基因組
建立索引：
cd /home/ngs/Pipeline/WES/database/gatk/hg38
gzip -d Homo_sapiens_assembly38.fasta.gz
mkdir index && cd index
nohup bwa index -a bwtsw -p hg38 ../Homo_sapiens_assembly38.fasta &
# -a有兩種構(gòu)建index的算法：bwtsw：大的基因組數(shù)據(jù)污朽，必須大于10MB崔梗，比如人的全基因組；is：默認(rèn)的算法柄慰，速度較快趟薄，需要較大的內(nèi)存贱枣，不能構(gòu)建大于2GB的數(shù)據(jù)庫
# -p str：輸出數(shù)據(jù)庫的前綴，默認(rèn)和輸入的文件名一致
進(jìn)行比對：
# 單樣本比對
INDEX=/home/ngs/Pipeline/WES/database/gatk/hg38/index/hg38
sample=025666
cd /home/ngs/Pipeline/WES/project/clean/
bwa mem -t 10 -R "@RG\tID:$sample\tSM:$sample\tLB:WGS\tPL:Illumina" $INDEX 025666_S168_L001.R1_val_1.fq.gz 025666_S168_L001.R2_val_2.fq.gz |samtools sort -@ 5 -o 025666.bam -
# -t：線程數(shù)
# -R接的是 Read Group的字符串信息日麸，它是用來將比對的read進(jìn)行分組的寄啼，這個信息對于后續(xù)對比對數(shù)據(jù)進(jìn)行錯誤率分析和Mark duplicate時非常重要逮光。ID是Read Group的分組ID，一般設(shè)置為測序的lane ID墩划；PL指的是所用的測序平臺涕刚；SM為樣本ID；LB為測序文庫的名字乙帮。這些信息設(shè)置好后杜漠，在RG字符串中要用制表符將它們分開
# 使用samtools將bwa產(chǎn)生的sam文件轉(zhuǎn)換為bam文件并排序（根據(jù)左起位點對序列排序），-@為線程數(shù)察净，后邊的“-”為管道輸出的占位符
# 批量比對
ls *R1_val_1.fq.gz >1
ls *R2_val_2.fq.gz >2
paste 1 2 > tmp
cat tmp |cut -f1|cut -d"_" -f1 >del
paste del tmp >config
INDEX=/home/ngs/Pipeline/WES/database/gatk/hg38/index/hg38
cat config |while read id
do
# 對config中的每一列進(jìn)行賦值驾茴，并且開始循環(huán)讀取和批量轉(zhuǎn)換
arr=($id)
fq1=${arr[1]} # fq1為第2列
fq2=${arr[2]} # fq2為第3列
sample=${arr[0]} #sample為第1列
bwa mem -t 10 -R "@RG\tID:$sample\tSM:$sample\tLB:WGS\tPL:Illumina" $INDEX $fq1 $fq2 |samtools sort -@ 5 -o $sample.bam -
done
