參考：
以水稻為例教你如何使用BSA方法進(jìn)行遺傳定位（上篇） - 簡(jiǎn)書(shū) (jianshu.com)
1. 數(shù)據(jù)質(zhì)控-fastp
#創(chuàng)建文件夾
mkdir 1.clean_data
cd clean_data
vi clean_data.sh
fastp -i ~/workspace/BSA/practice/data/SRR6327815_1.fastq.gz -I ~/workspace/BSA/practice/data/SRR6327815_2.fastq.gz -o SRR6327815_1.fastq.gz -O SRR6327815_2.fastq.gz
fastp -i ~/workspace/BSA/practice/data/SRR6327816_1.fastq.gz -I ~/workspace/BSA/practice/data/SRR6327816_2.fastq.gz -o SRR6327816_1.fastq.gz -O SRR6327816_2.fastq.gz
fastp -i ~/workspace/BSA/practice/data/SRR6327817_1.fastq.gz -I ~/workspace/BSA/practice/data/SRR6327817_2.fastq.gz -o SRR6327817_1.fastq.gz -O SRR6327817_2.fastq.gz
fastp -i ~/workspace/BSA/practice/data/SRR6327818_1.fastq.gz -I ~/workspace/BSA/practice/data/SRR6327818_2.fastq.gz -o SRR6327818_1.fastq.gz -O SRR6327818_2.fastq.gz
chmod +x clean_data.sh
sh clean_data.sh &
2. 序列比對(duì)
①參考文章中的shell腳本
# BSA/practice/目錄下
mkdir -p 2.ref_align
#給變量賦值
NUMBER_THREADS=12
REFERENCE=data/ref/IRGSP-1.0_genome.fasta
# for循環(huán)
for i in `ls 1.clean_data/`; do
#{i%%_*}:刪去i的第一個(gè)"_"及其右邊字符串
sample=${i%%_*};
(bwa mem -M -R "@RG\\tID:${sample}\\tSM:${sample}\\tPL:ILLUMINA" -t $NUMBER_THREADS $REFERENCE 1.clean_data/${sample}_1.fastq.gz 1.clean_data/${sample}_2.fastq.gz || echo -n 'error' )
#排序
| samtools sort -@ 20 -o 2.ref_align/${sample}_sort.bam -;
#構(gòu)建索引
samtools index -@ 20 2.ref_align/${sample}_sort.bam;
done
②修改版
參考文章里的這一段命令其實(shí)相當(dāng)于bwa比對(duì)同一個(gè)文件比對(duì)兩次
#可以看輸出的sample 更直觀點(diǎn)
for i in `ls 1.clean_data/`;do sample=${i%%_*}; echo $sample;done
image.png
所以我修改了一下~
#先將這個(gè)sample重定向輸出到一個(gè)文件中
for i in `ls 1.clean_data/`;do
sample=${i%%_*};
echo $sample >> sample.txt;
done
uniq sample.txt>sample.txt
cat sample.txt|while read line
do
bwa mem -t 16 -M -Y -R "@RG\tID:$line\tPL:ILLUMINA\tLB:$line\tSM:$line" $REFERENCE ../1.clean_data/${line}_1.fastq.gz ../1.clean_data/${line}_2.fastq.gz >${line}.bam &&
# 排序
samtools sort -@ 16 -m 4G -O bam ./${line}.bam -o ./${line}_sort.bam && \
samtools index ./${line}_sort.bam
done
3. 去除重復(fù)序列
去除相同位置且序列一模一樣的reads，通常是由于PCR擴(kuò)增引起。PCR擴(kuò)增中也有可能會(huì)出現(xiàn)錯(cuò)誤习柠，序列擴(kuò)增錯(cuò)誤可能會(huì)被變異檢測(cè)軟件判定為變異蛉艾，過(guò)濾掉重復(fù)序列有利于提升后續(xù)檢測(cè)的準(zhǔn)確性。
sambamba速度比較快晤愧，且結(jié)果和picard一模一樣大莫。
cat sample.txt|while read line
do
sambamba markdup -r -t 10 ${line}_sort.bam ${line}_mkdup.bam
done
查看bam文件可使用如下命令
samtools view -h SRR6327815_mkdup.bam|less -S
4. 變異檢測(cè)
變異檢測(cè)常用軟件有bcftools, freebayes和GATK.
參考文章里用到的是bcftools，主要原因還是它的速度比較快官份。
為了讓他的速度更快只厘，使用--region參數(shù)分染色體并行，由于水稻有12條染色體舅巷，相當(dāng)于提速了12倍
# BSA/practice目錄下
mkdir -p 3.variants
cd 3.variants
for i in `ls ../2.ref_align/*_mkdup.bam`;do
sample=${i##*/};#去掉最后一個(gè)/及其左邊的字符串
echo $sample >> bam_list.txt;
done
#seqkit seq -n 羔味，name，輸出染色體名稱(chēng)
seqkit seq -n ../data/ref/IRGSP-1.0_genome.fasta | \
while read region
do
bcftools mpileup -f ../data/ref/IRGSP-1.0_genome.fasta \
--redo-BAQ --min-BQ 30 \
--per-sample-mF \
--annotate FORMAT/AD,FORMAT/DP \
--regions ${region} \
-Ou --bam-list bam_list.txt | \
bcftools call -mv -Ob -o 03-variants/${region}.bcf &
done
接著將拆分結(jié)果合并
# BSA/practice目錄下
mkdir -p 4.variants_filter
cd 4.variants_filter
bcftools concat --naive -o merged.bcf ../3.variants/*.bcf
5. 變異過(guò)濾
根據(jù)需求進(jìn)行不同閾值或條件設(shè)置來(lái)對(duì)變異檢測(cè)結(jié)果進(jìn)行過(guò)濾
一般需要考慮的因素有：
①測(cè)序深度钠右；
②非參考基因組的高質(zhì)量reads數(shù)介评；
③是否和indel緊鄰，通常和indel比較近的snp都不可靠爬舰；
④平均的比對(duì)質(zhì)量们陆。
參考文章中是這樣
bcftools filter -g3 -G10 -e'%QUAL<10 || (RPB<0.1 && %QUAL<15) || (AC<2 && %QUAL<15) || MQ < 30 || MQSB <=0.1' merged.bcf > filter.vcf
但是不知道為什么報(bào)錯(cuò)如下
[filter.c:2903 filters_init1] Error: the tag "RPB" is not defined in the VCF header
[filter.c:2903 filters_init1] Error: the tag "MQSB" is not defined in the VCF header
所以我把RPB和MQSB的過(guò)濾條件給刪掉了，最后用的如下命令
bcftools filter -g3 -G10 -e'%QUAL<10 || (AC<2 && %QUAL<15) || MQ < 30 ' merged.bcf > filter.vcf
PS：后來(lái)我查看了一下文件情屹，里面是RPBZ和MQSBZ
image.png
bcftools filter -g3 -G10 -e'%QUAL<10 || (RPBZ<0.1 && %QUAL<15) || (AC<2 && %QUAL<15) || MQ < 30 || MQSBZ <=0.1' merged.bcf > full_filter.vcf
#對(duì)比一下完整過(guò)濾條件和部分過(guò)濾條件后的SNP數(shù)
grep -v "#" filter.vcf| wc -l
1037996
grep -v "#" full_filter.vcf| wc -l
266877
但是這最后與參考文章找那個(gè)差的也有點(diǎn)太多了坪仇。。
接著選擇SNP來(lái)進(jìn)行后續(xù)的分析
bcftools view -i 'TYPE="snp" & N_ALT =1 & STRLEN(ALT) = 1' filter.vcf > snps.vcf
#統(tǒng)計(jì)一下有多少SNPs
grep -v "#" snps.vcf |wc -l
#顯示有90w垃你，參考文章中是80w椅文，可能是過(guò)濾條件那里我刪去了兩個(gè)
#我決定先繼續(xù)后面的數(shù)據(jù)分析
904821
