03、測(cè)序數(shù)據(jù)批量比對(duì)到參考基因組

建立索引：

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ù)庫(kù)
# -p str：輸出數(shù)據(jù)庫(kù)的前綴，默認(rèn)和輸入的文件名一致

進(jìn)行比對(duì)：

# 單樣本比對(duì)
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的字符串信息，它是用來(lái)將比對(duì)的read進(jìn)行分組的，這個(gè)信息對(duì)于后續(xù)對(duì)比對(duì)數(shù)據(jù)進(jìn)行錯(cuò)誤率分析和Mark duplicate時(shí)非常重要。ID是Read Group的分組ID，一般設(shè)置為測(cè)序的lane ID；PL指的是所用的測(cè)序平臺(tái)；SM為樣本ID；LB為測(cè)序文庫(kù)的名字。這些信息設(shè)置好后，在RG字符串中要用制表符將它們分開
# 使用samtools將bwa產(chǎn)生的sam文件轉(zhuǎn)換為bam文件并排序（根據(jù)左起位點(diǎn)對(duì)序列排序），-@為線程數(shù)，后邊的“-”為管道輸出的占位符

# 批量比對(duì)
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
# 對(duì)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