注：本次使用的数据为不同品种李子重测序
1.准备工作
首先下载相应的参考基因组及注释文件{本次使用的参考基因组为三月李的，参考基因组从蔷薇科植物基因组数据库下载（GDR)}
安装好GATK,samtools,fastp,bwa,vcftools,以上软件均通过conda安装
参考基因组存放位置
重测序原始数据存放位置
1.1 使用samtools构建参考基因组索引
samtools faidx referencegenome.fa #索引文件后缀.fai
1.2 使用bwa构建参考基因组索引
bwa index referencegenome.fa #bwa索引文件有.amb,ann,bwt,pac,sa
1.3 使用gatk构建参考基因组dict文件
time gatk CreateSequenceDictionary -R referencegenome.fa -O referencegenome.dict
-R 参考基因组绝对路径 -O 输出文件路径
2.原始数据质控、比对、排序标记PCR重复
2.1原始数据质控
ls *.fq.gz|paste - - |perl -lane 'print qq{fastp -w 4 -i $F[0] -I $F[1] -o $F[0].clean -O $F[1].clean}'|perl -lpe 's/.fq.gz.clean/.clean.fq.gz/g' | while read line;do $line;done
2.2 比对和排序
2.2.1 创建将待处理文件名整合为一个文件
ls *clean.fq.gz >pairend1
ls *clean.fq.gz >pairend2
paste pairend1 pairend2 >bwa.config
将第一列设为样本名并以tab键分隔
2.2.2 创建脚本文件
vim bwa.sh 文件内容如下
#!/bin/bash
#SBATCH -p amd_512 ##指定使用的队列名称
#SBATCH -N 1 ##指定使用的节点数是1
#SBATCH -n 1 ##指定使用的核数为1（串行的话使用1核，如果需要使用并行可以修改核数）
#SBATCH --mem-per-cpu=8G #后台可能因为cpu不够而将程序杀掉可将cpu设置的相对大点
sample="DBL FTL HFTL QPFTL ML1 SIYL WSL WSQCL" # 样品名根据自己项目更改
echo $sample
INDEX=/public1/home/scg2749/project/Pruns/-----/X101SC22043754-Z01-J002/04.Reference/Prunus_salicina_Sanyueli_FAAS_v1.0.fa #参考基因组位置（你的参考基因组存放位置，最好是绝对路径）
cat bwa.config |while read id
do
arr=($id)
fq1=${arr[1]}
fq2=${arr[2]}
sample=${arr[0]}
bwa mem -K 100000000 -Y -t 8 -R "@RG\tID:$sample\tSM:$sample\tLB:WGS\tPL:Illumina" $INDEX $fq1 $fq2|
samtools sort -@ 16 -o /public1/home/scg2749/project/Pruns/xieweiwei/X101SC22043754-Z01-J002/02.bam/$sample.bam -
done
提交任务 sbatch bwa.sh 查看任务squeue 取消任务 scancel job ID
2.3 标记PCR重复
此时使用的是GATK子程序MarkDuplicates
创建MarkDuplicates.sh内容如下
#!/bin/bash
#SBATCH -p amd_512 ##指定使用的队列名称
#SBATCH -N 4 ##指定使用的节点数是1
#SBATCH -n 4 ##指定使用的核数为1（串行的话使用1核，如果需要使用并行可以修改核数）
#SBATCH --mem-per-cpu=16G
ls *.bam|perl -lane 'print qq{time gatk MarkDuplicates -I /public1/home/scg2749/project/Pruns/----/X101SC22043754-Z01-J002/02.bam/$F[0] -O /public1/home/scg2749/project/Pruns/----/X101SC22043754-Z01-J002/03.sorted.bam/$F[0].markdup.bam -M $F[0].sorted.markdup_metrics.txt}'|perl -lpe 's/.bam.markdup.bam/.markdup.bam/'|perl -lpe 's/.bam.sorted.markdup_metrics.txt/.markdup_metrics.txt/'|while read line;do $line;done
提交任务 sbatch MarkDuplicates.sh如果不是在服务器上科sh MarkDuplicates.sh提交任务
2.4 对上步得到的bam文件创建索引
sbatch samtools_index.sh
#!/bin/bash
#SBATCH -p amd_512 ##指定使用的队列名称
#SBATCH -N 1 ##指定使用的节点数是1
#SBATCH -n 1 ##指定使用的核数为1（串行的话使用1核，如果需要使用并行可以修改核数）
#SBATCH --mem-per-cpu=16G
ls *.bam|perl -lane 'print qq{time samtools index /public1/home/scg2749/project/Pruns/xieweiwei/X101SC22043754-Z01-J002/03.sorted.bam/$F[0]}'|while read line;do $line;done
3.使用GATK检测变异位点
3.1 先将bam文件生成中间文件g.vcf
此时使用gatk HaplotypeCaller子程序
vim HaplotypeCaller.sh
#!/bin/bash
#SBATCH -p amd_512 ##指定使用的队列名称
#SBATCH -N 4 ##指定使用的节点数是1
#SBATCH -n 4 ##指定使用的核数为1（串行的话使用1核，如果需要使用并行可以修改核数）
#SBATCH --mem-per-cpu=8G
ls *.bam|perl -lane 'print qq{time gatk HaplotypeCaller -R /public1/home/scg2749/project/Pruns/----/X101SC22043754-Z01-J002/04.Reference/Prunus_salicina_Sanyueli_FAAS_v1.0.fa --emit-ref-confidence GVCF -I /public1/home/scg2749/project/Pruns/----/X101SC22043754-Z01-J002/03.sorted.bam/$F[0] -O /public1/home/scg2749/project/Pruns/----/X101SC22043754-Z01-J002/05.result/$F[0].g.vcf}'|perl -lpe 's/.rmdup.bam.g.vcf/.g.vcf/'|while read line;do $line;done
sbatch HaplotypeCaller.sh
参数详解 -R 参考基因组位置 -I bam文件位置 -O输出路径 注意大小写
3.2使用GenotypeGVCFs对GVCF进行joint calling 通过gvcf检测变异
vim GenotypeGVCFs.sh
#!/bin/bash
#SBATCH -p amd_512 ##指定使用的队列名称
#SBATCH -N 4 ##指定使用的节点数是1
#SBATCH -n 4 ##指定使用的核数为1（串行的话使用1核，如果需要使用并行可以修改核数)
ls *.g.vcf|perl -lane 'print qq{time gatk GenotypeGVCFs -R /public1/home/scg2749/project/Pruns/xieweiwei/X101SC22043754-Z01-J002/04.Reference/Prunus_salicina_Sanyueli_FAAS_v1.0.fa -V /public1/home/scg2749/project/Pruns/xieweiwei/X101SC22043754-Z01-J002/05.result/$F[0] -O /public1/home/scg2749/project/Pruns/xieweiwei/X101SC22043754-Z01-J002/05.result/$F[0].vcf}'|perl -lpe 's/.markdup.bam.g.vcf.vcf/.vcf/'g|while read line;do $line;done
sbatch GenotypeGVCFs.sh
参数详解：-R 参考基因组位置 -V 上一步中生成的g.vcf文件路径 -O输出路径
为节省储存可使用gzip对vcf文件压缩
ls *.vcf|perl -lane 'print qq{time bgzip -f /public1/home/scg2749/project/Pruns/----/X101SC22043754-Z01-J002/05.result/$F[0]}'|while read line;do $line;done
构建tabix索引
ls *.vcf.gz|while read id;do time tabix -p vcf /public1/home/scg2749/project/Pruns/xieweiwei/X101SC22043754-Z01-J002/05.result/$id;done
此时就得到不同样本相对于参考基因组的变异情况了，后序还需对变异数据进行质控及注释
写在最后，
一天过去了，又划了一天水。最近学校食堂开了新窗口了是瓦罐汤和拌粉，开业第一天去尝了一下，从老板的动作来看属于是第一次开店，但是手艺还不错至少我觉得拌粉挺好吃的，希望老板能长久开下去吧。感觉自己最近失去的挺多的，包括人和事，心情挺糟糕的。打算看看文献看能不能分散注意力，顺带做做文献分享，看看自己行不行，毕竟湿实验也做不出来。
网友评论