步骤1：数据质控：

nohup fastqc -t 24 * &
mkdir filter

fastp代码如下，保存为fastp.sh

#!/bin/bash

for file in `ls *_R1.fq.gz | perl -lpe 's/_R1.fq.gz//'`; do 
    fastp -i ${file}_R1.fq.gz -f 5 -t 0 -o ./filter/${file}_fp_R1.fq.gz \
    -I  ${file}_R2.fq.gz -F 5 -T 30  -O ./filter/${file}_fp_R2.fq.gz \
    -h  ${file}_report.html \
    -j  ${file}_report.json 

done

执行fastp命令，批量处理fastq文件

sh fastp.sh

步骤2：call SNP，使用GATK的流程
genome文件构建索引

bwa index public.fasta
samtools  faidx public.fasta
gatk CreateSequenceDictionary -R public.fasta -O public.dict

生成每个样本的gvcf文件，限制了内存和线程数目

#! /bin/sh



for sample in `cat accession1.txt`;do
       bwa mem -t 8 -M -R "@RG\tID:${sample}\tSM:${sample}\tLB:WGS\tPL:Illumina" ../genome/NAU.genome.fa ../data/${sample}_fp_R1.fq.gz ../data/${sample}_fp_R2.fq.gz | samtools view -Sb - > ${sample}.bam && echo "**bwa mapping done **"
       samtools sort -@ 8 -o ${sample}.sorted.bam ${sample}.bam
       rm ${sample}.bam

        #mark duplicate
        gatk --java-options "-Xmx48G -XX:ParallelGCThreads=8 -Djava.io.tmpdir=./tmp" MarkDuplicates -I ${sample}.sorted.bam -O ${sample}.sorted.markdup.bam -M ${sample}.sorted.markdup_metrics.txt  --ASSUME_SORTED true
        

        
        #index
        samtools index ${sample}.sorted.markdup.bam

        #statistic
        #samtoolsflagstat -@ 24 ${sample}.sorted.markdup.bam >${sample}.sorted.markdup.stat
#       qualimap --java-mem-size=30G  bamqc -bam ${sample}.sorted.markdup.bam

        #call variant
        gatk --java-options "-Xmx48G"  HaplotypeCaller --native-pair-hmm-threads 8  -R ../genome/NAU.genome.fa -I ${sample}.sorted.markdup.bam -ERC GVCF -O ${sample}.erc.gvcf

        rm ${sample}.sorted.bam
        rm ${sample}.sorted.markdup.bam
done

合并gvcf

#! /bin/sh
gatk --java-options "-Xmx48G" CombineGVCFs -R ../genome/NAU.genome.fa $(for i in $(ls *.gvcf);do echo "--variant $i";done) -O combined.gvcf
gatk --java-options "-Xmx48G" GenotypeGVCFs -R ../genome/NAU.genome.fa -V combined.gvcf -O combined.vcf

SNP 过滤

pop=$1
## reference genome
reference=$2

### software source
## GATK: https://gatk.broadinstitute.org/hc/en-us

### Steps
## select SNPs from raw vcf
gatk SelectVariants --select-type-to-include SNP --reference ${reference} --variant ${pop}.vcf.gz --output ${pop}.raw.snp.vcf.gz

## perform hard filtering
gatk VariantFiltration --variant ${pop}.raw.snp.vcf.gz --filter-expression "QD < 2.0 || FS > 60.0 || MQ < 40.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "snp_filter" --genotype-filter-expression "DP < 2 || DP > 50" --genotype-filter-name "dp_fail" --output ${pop}.flag.snp.vcf.gz

## selecting bi-allelic SNPs that pass the filtering
gatk SelectVariants --exclude-filtered true --restrict-alleles-to BIALLELIC --reference ${reference} --variant ${pop}.flag.snp.vcf.gz --output ${pop}.hardfiltered.snp.vcf.gz

## perform population filtering
vcftools --gzvcf ${pop}.hardfiltered.snp.vcf.gz --max-missing 0.9 --maf 0.05 --recode --recode-INFO-all --stdout | bgzip -c >  ${pop}.maffiltered.snp.vcf.gz

QTL-seq，参考https://laowang2023.cn/2023/11/07/20231107-easyQTLseq/
使用 GATK VariantsToTable会报错，改成以下命令，gatk版本是4.0

gatk IndexFeatureFile -I /home/data/t010537/wjl-BSA/gvcf/combined.maffiltered.snp.vcf.gz
gatk  --java-options "-Xmx30g" VariantsToTable -F CHROM -F POS -F REF -F ALT -GF GT -GF AD -GF GQ  -V combined.maffiltered.snp.vcf.gz -O BSA.filter.SNPs.table

使用easyQTL

library(easyQTLseq)
file_path <- "BSA.filter.SNPs.table"
data <- readr::read_tsv(file = file_path)
data <- data[!grepl("^sca", data$CHROM), ]
x <- select_sample_and_SNP(data = data, highP = "D60062-343_fp", lowP = "D60062-344_fp", highB = "D60062-chunchi3_fp", lowB = "D60062-chunchi4_fp", popType = "F2", bulkSize = c(30, 30))
x_filter <- filterDP(x = x)
x_filter <- calc_index_etc(x = x_filter, outPrefix = "outprefix", winSize = 2000000, winStep = 20000)

export_figure(x = x_filter, 
              outPrefix = "outprefix", 
              minN = 20,              
              width = 6, height = 3)