1. 数据质控-fastp
#创建文件夹
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. 序列比对
①参考文章中的shell脚本
# BSA/practice/目录下
mkdir -p 2.ref_align
#给变量赋值
NUMBER_THREADS=12
REFERENCE=data/ref/IRGSP-1.0_genome.fasta
# for循环
for i in `ls 1.clean_data/`; do
#{i%%_*}:删去i的第一个"_"及其右边字符串
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 -;
#构建索引
samtools index -@ 20 2.ref_align/${sample}_sort.bam;
done
②修改版
参考文章里的这一段命令其实相当于bwa比对同一个文件比对两次
#可以看输出的sample 更直观点
for i in `ls 1.clean_data/`;do sample=${i%%_*}; echo $sample;done
image.png
所以我修改了一下~
#先将这个sample重定向输出到一个文件中
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. 去除重复序列
去除相同位置且序列一模一样的reads，通常是由于PCR扩增引起。PCR扩增中也有可能会出现错误，序列扩增错误可能会被变异检测软件判定为变异，过滤掉重复序列有利于提升后续检测的准确性。
sambamba速度比较快，且结果和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. 变异检测
变异检测常用软件有bcftools, freebayes和GATK.
参考文章里用到的是bcftools，主要原因还是它的速度比较快。
为了让他的速度更快，使用--region参数分染色体并行，由于水稻有12条染色体，相当于提速了12倍
# BSA/practice目录下
mkdir -p 3.variants
cd 3.variants
for i in `ls ../2.ref_align/*_mkdup.bam`;do
sample=${i##*/};#去掉最后一个/及其左边的字符串
echo $sample >> bam_list.txt;
done
#seqkit seq -n ，name，输出染色体名称
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
接着将拆分结果合并
# BSA/practice目录下
mkdir -p 4.variants_filter
cd 4.variants_filter
bcftools concat --naive -o merged.bcf ../3.variants/*.bcf
5. 变异过滤
根据需求进行不同阈值或条件设置来对变异检测结果进行过滤
一般需要考虑的因素有：
①测序深度；
②非参考基因组的高质量reads数；
③是否和indel紧邻，通常和indel比较近的snp都不可靠；
④平均的比对质量。
参考文章中是这样
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
但是不知道为什么报错如下
[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的过滤条件给删掉了，最后用的如下命令
bcftools filter -g3 -G10 -e'%QUAL<10 || (AC<2 && %QUAL<15) || MQ < 30 ' merged.bcf > filter.vcf
PS：后来我查看了一下文件，里面是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
#对比一下完整过滤条件和部分过滤条件后的SNP数
grep -v "#" filter.vcf| wc -l
1037996
grep -v "#" full_filter.vcf| wc -l
266877
但是这最后与参考文章找那个差的也有点太多了。。
接着选择SNP来进行后续的分析
bcftools view -i 'TYPE="snp" & N_ALT =1 & STRLEN(ALT) = 1' filter.vcf > snps.vcf
#统计一下有多少SNPs
grep -v "#" snps.vcf |wc -l
#显示有90w，参考文章中是80w，可能是过滤条件那里我删去了两个
#我决定先继续后面的数据分析
904821
