墨滴

生信探索

V1

2023/03/31阅读：28主题：姹紫

基因组实战03: WGS toy example

<生信交流与合作请关注公众号@生信探索>

借鉴Reference中第2、3篇文章的代码。分析的数据是大肠杆菌，因为基因组小，适合拿来快速跑通整个流程

00 下载fastq数据

mkdir -p ~/Project/DNA/raw
cd ~/Project/DNA/raw
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR177/003/SRR1770413/SRR1770413_1.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR177/003/SRR1770413/SRR1770413_2.fastq.gz

01 filter the bad quality reads and remove adaptors

  • 脚本 
#!/bin/bash
#clean.sh
# python3.10/python2.7
# fastqc multiqc trim_galore

HOME_DIR=~/Project/DNA
FASTQ_DIR=${HOME_DIR}/raw
CLEAN_DIR=${HOME_DIR}/clean
N_JOBS=12

if  [ ! -d $CLEAN_DIR ]
then
    mkdir -p $CLEAN_DIR
fi

micromamba activate dna2

# 1.QC
cd $FASTQ_DIR
fastqc --thread $N_JOBS *fastq.gz && multiqc .

# 2.remove adapter
micromamba activate dna3
for FILE in `ls $FASTQ_DIR/*_1.fastq.gz`
do
    SAMPLE=`basename $FILE | sed s/_1\.fastq\.gz//`
    echo "Processing $SAMPLE"

    F1=$FASTQ_DIR/$SAMPLE"_1.fastq.gz"
    F2=$FASTQ_DIR/$SAMPLE"_2.fastq.gz"

    trim_galore \
        --quality 30 \
        --length 50 \
        --output_dir $CLEAN_DIR \
        --cores $N_JOBS \
        --paired \
        --fastqc \
        -e 0.1 \
        --trim-n \
        $F1 $F2

    mv $CLEAN_DIR/${SAMPLE}"_1_val_1.fq.gz" \
        $CLEAN_DIR/${SAMPLE}"_1.fastq.gz"

    mv $CLEAN_DIR/${SAMPLE}"_2_val_2.fq.gz" \
        $CLEAN_DIR/${SAMPLE}"_2.fastq.gz"

    mv $CLEAN_DIR/${SAMPLE}"_1_val_1_fastqc.html" \
        $CLEAN_DIR/${SAMPLE}"_1_fastqc.html"

    mv $CLEAN_DIR/${SAMPLE}"_1_val_1_fastqc.zip" \
        $CLEAN_DIR/${SAMPLE}"_1_fastqc.zip"

    mv $CLEAN_DIR/${SAMPLE}"_2_val_2_fastqc.html" \
        $CLEAN_DIR/${SAMPLE}"_2_fastqc.html"

    mv $CLEAN_DIR/${SAMPLE}"_2_val_2_fastqc.zip" \
        $CLEAN_DIR/${SAMPLE}"_2_fastqc.zip"
done

# 3.QC for clean
cd $CLEAN_DIR
micromamba activate dna2
multiqc .

  • 运行 
source clean.sh &> clean.sh.log &

02 align

下载参考基因组数据

curl -OJX GET "https://api.ncbi.nlm.nih.gov/datasets/v2alpha/genome/accession/GCF_000005845.2/download?include_annotation_type=GENOME_FASTA,GENOME_GFF,RNA_FASTA,CDS_FASTA,PROT_FASTA,SEQUENCE_REPORT&filename=GCF_000005845.2.zip" -H "Accept: application/zip"
unzip GCF_000005845.2.zip
cp ncbi_dataset/data/GCF_000005845.2/GCF_000005845.2_ASM584v2_genomic.fna genome/E.coli.fa

生成的文件索引

#!/bin/bash
#python2.7
HOME_DIR=~/Project/DNA
GENOME_FILE=E.coli.fa
INDEX_DIR=$HOME_DIR/genome
micromamba activate dna2
echo "----------------- BWA INDEX START-----------------"
bwa index $INDEX_DIR/$GENOME_FILE && \
echo "----------------- BWA INDEX DONE -----------------"

echo "----------------- CreateSequenceDictionary START -----------------"
gatk CreateSequenceDictionary \
    --REFERENCE $INDEX_DIR/$GENOME_FILE \
    --OUTPUT $INDEX_DIR/E.coli.dict && echo "** dict done **"  && \
echo "----------------- CreateSequenceDictionary DONE -----------------"

samtools faidx $INDEX_DIR/$GENOME_FILE
source index.sh &> index.sh.log &

对比

-R：Read group；

ID：Read group identifier，必须唯一，一般为样本名；

SM：Sample，样本名；

LB：Library，测序文库；

PL：Platform，测序平台

#!/bin/bash
#python2.7
micromamba activate dna2
HOME_DIR=~/Project/DNA
BWA_INDEX=$HOME_DIR/genome/E.coli.fa
CLEAN_DIR=$HOME_DIR/clean
ALIGN_DIR=$HOME_DIR/align
N_JOBS=12

if  [ ! -d $ALIGN_DIR ]
then
    mkdir -p $ALIGN_DIR
fi

for FILE in `ls $CLEAN_DIR/*_1.fastq.gz`
do
    SAMPLE=`basename $FILE | sed s/_1\.fastq\.gz//`

    echo "----------------- BWA ${SAMPLE} START -----------------"
    F1=$CLEAN_DIR/$SAMPLE"_1.fastq.gz"
    F2=$CLEAN_DIR/$SAMPLE"_2.fastq.gz"
    RG="@RG\tID:"$SAMPLE"\tSM:"$SAMPLE"\tLB:WGS\tPL:ILLUMINA"
    bwa mem -t $N_JOBS -R $RG $BWA_INDEX $F1 $F2 | \
        samtools sort --threads $N_JOBS -O bam -o $ALIGN_DIR/$SAMPLE".bam" - && \
    echo  "----------------- BWA ${SAMPLE} DONE -----------------"

    echo "----------------- MarkDuplicates ${SAMPLE} START-----------------"
    gatk MarkDuplicates \
        --INPUT ${ALIGN_DIR}/${SAMPLE}".bam" \
        --OUTPUT ${ALIGN_DIR}/${SAMPLE}".mark.bam" \
        --METRICS_FILE ${ALIGN_DIR}/${SAMPLE}".mark.matrics" && \
    echo  "----------------- MarkDuplicates ${SAMPLE} DONE -----------------" &&
    rm -f ${ALIGN_DIR}/${SAMPLE}".bam"

    echo "-----------------Samtools Indexing ${SAMPLE} START -----------------"
    samtools index ${ALIGN_DIR}/${SAMPLE}".mark.bam"  && \
    echo  "----------------- Samtools Indexing ${SAMPLE} DONE -----------------"
done

source bwa.sh &> bwa.sh.log &

03 call variants

#!/bin/bash
micromamba activate dna2
HOME_DIR=~/Project/DNA
BWA_INDEX=$HOME_DIR/genome/E.coli.fa
ALIGN_DIR=$HOME_DIR/align
MUTATION=$HOME_DIR/mutation
N_JOBS=12

if  [ ! -d $MUTATION ]
then
    mkdir -p $MUTATION
fi

# 1 每个样本生成一个gvcf文件
for FILE in `ls ${ALIGN_DIR}/*.mark.bam`
do
    SAMPLE=`basename $FILE | sed s/\.mark\.bam//`
    echo $SAMPLE

    gatk HaplotypeCaller \
    --reference $BWA_INDEX \
    --emit-ref-confidence GVCF \
    --input $ALIGN_DIR/$SAMPLE".mark.bam" \
    --output $MUTATION/$SAMPLE".gvcf" && echo "** gvcf done **"
done

# 2 通过gvcf检测变异
gatk GenotypeGVCFs \
--reference $BWA_INDEX \
--variant $MUTATION/$SAMPLE".gvcf" \
--output $MUTATION/final.vcf && echo "** vcf done **"

# 3 压缩 构建tabix索引
bgzip -f $MUTATION/final.vcf
tabix -p vcf $MUTATION/final.vcf.gz

# 4.过滤（硬指标）

## 使用SelectVariants，选出SNP
gatk SelectVariants \
    -select-type SNP \
    --variant $MUTATION/final.vcf.gz \
    --output $MUTATION/final.snp.vcf.gz

## 为SNP作硬过滤
gatk VariantFiltration \
    --variant $MUTATION/final.snp.vcf.gz \
    --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
    --filter-name "PASS" \
    --output $MUTATION/final.snp.filter.vcf.gz

## 使用SelectVariants，选出Indel
gatk SelectVariants \
    -select-type INDEL \
    --variant $MUTATION/final.vcf.gz \
    --output $MUTATION/final.indel.vcf.gz

## 为Indel作过滤
gatk VariantFiltration \
    --variant $MUTATION/final.indel.vcf.gz \
    --filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
    --filter-name "PASS" \
    --output $MUTATION/final.indel.filter.vcf.gz

# 重新合并过滤后的SNP和Indel
gatk MergeVcfs \
    --INPUT $MUTATION/final.snp.filter.vcf.gz \
    --INPUT $MUTATION/final.indel.filter.vcf.gz \
    --OUTPUT $MUTATION/final.filter.vcf.gz

cd $MUTATION
rm -f *snp* *indel* final.vcf* *gvcf*
source call_variant.sh &> call_variant.sh.log &

Reference

# gatk api
https://gatk.broadinstitute.org/hc/en-us/articles/13832655155099--Tool-Documentation-Index
#GATK4.0和全基因组数据分析实践  
https://mp.weixin.qq.com/s/Heu-jmJeM3EQy2PmyA-gyQ
https://mp.weixin.qq.com/s/ooU7Tuh0adGNKEISELFjUA
#GATK best practices
https://mp.weixin.qq.com/s/GvuL8NlFSkQ6MVdghlEYUQ
# trim-galore
https://mp.weixin.qq.com/s/n3G_qot5mhB3ZR1gcDnV6g
# 方法分享 | 全外显子测序分析的标准流程
https://mp.weixin.qq.com/s/WOeLHXIomhNjSFXjXZ8zcg
https://mp.weixin.qq.com/s/WOeLHXIomhNjSFXjXZ8zcg
# GATK4全基因组数据分析最佳实践  
https://mp.weixin.qq.com/s/r3sghKtEKq1FnjQ05WCcnA
# 第1节 测序技术  
https://mp.weixin.qq.com/s/kq2i5tNZmm9GKGAx2Day-g
#第2节 FASTA和FASTQ  
https://mp.weixin.qq.com/s/0h5B3T6RKTHTsZBuoZvtgQ
#第3节 数据质控  
https://mp.weixin.qq.com/s/8hY0U48kiVTH6Yx4JBcAhg#
#第4节 构建WGS主流程  
https://mp.weixin.qq.com/s/AT2oodvdqWgbj1gvVo1hWQ
#第5节 理解并操作BAM文件  
https://mp.weixin.qq.com/s/UnMyAuUHmK7DMGo8h88oQw

# WGS(全基因组测序)/WES(全外显子组测序)/WGRS(全基因组重测序)分析与可视化教程 
https://zhuanlan.zhihu.com/p/380684876

https://blog.csdn.net/bio_meimei/article/details/108236406

https://zhuanlan.zhihu.com/p/422365954

https://zhuanlan.zhihu.com/p/69726572
https://hpc.nih.gov/training/gatk_tutorial/haplotype-caller.html

分类：

后端

标签：

大数据

作者介绍

生信探索
V1

微信公众号：生信探索