WGS分析实战-01：从SRA数据下载到构建GenomicsDatabase

# mkdir 00.raw_data && mkdir 00.ref
prefetch --option-file SraAccList.txt

# cd ../00.ref
wget http://ftp.ensemblgenomes.org/pub/plants/release-52/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
gunzip Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz --stdout > arab_ref.fa

samtools faidx arab_ref.fa
# 结果文件：arab_ref.fa.fai
bwa index arab_ref.fa
# 结果文件：
# arab_ref.fa.amb  
# arab_ref.fa.ann  
# arab_ref.fa.bwt   
# arab_ref.fa.pac  
# arab_ref.fa.sa

# 工作目录：00.raw_data
# 存放到一个文件夹下
mkdir sra_file
for id in `ls -d D*`;do mv ${id}/${id}.sra sra_file;done

# 格式转换
ls sra_file/* > sra_id.txt
for id in `cat sra_id.txt`;do echo "fastq-dump --split-3 --gzip -O ./ ${id}" >> fastq-dump.sh;done

# ParaFly运行
ParaFly -c fastq-dump.sh -CPU 16 2>fastqdump.err.log
# rm -rf sra_file/

# 工作目录；00.raw_data
# mkdir ../01.clean_data
cat SraAccList.txt | while read id
do
 echo "fastp -i ${id}_1.fastq.gz -I ${id}_2.fastq.gz -o ../01.clean_data/${id}_clean_1.fastq.gz -O ../01.clean_data/${id}_clean_2.fastq.gz --cut_tail --n_base_limit 3 --length_required 60 --correction --thread 4 --html ../01.clean_data/${id}.html --json ../01.clean_data/${id}.json" >> fastp.commandlines
done

# ParaFly运行
ParaFly -c fastp.commandlines -CPU 16 2>fastp.err.log &

# mkdir ../02.bwa
cat ../00.raw_data/SraAccList.txt | while read id;do echo "bwa mem -t 6 -R '@RG\tID:foo\tPL:illumina\tSM:${id}' ../00.ref/arab_ref.fa ${id}_clean_1.fastq.gz ${id}_clean_2.fastq.gz | samtools view -Sb - > ../02.bwa/${id}.unsorted.bam" >> bwa.commondlines; done
# ParaFly运行
ParaFly -c bwa.commondlines -CPU 8 2>bwa_mem.err.log &

ls *.bam > unsorted_bam_id.txt
cat unsorted_bam_id.txt | while read id
do
 name=${id%%.*}
 echo "samtools sort -@ 6 -m 4G -O bam -o ${name}.sorted.bam ${id}" >> sortbam.commandlines
done

# ParaFly运行
ParaFly -c sortbam.commandlines -CPU 16 2>sortbam.err.log &
# 删除unsorted
# ls * | grep "unsorted.bam$" | xargs -i rm -rf {}

ls *.sorted.bam > sorted_bam_id.txt
cat sorted_bam_id.txt | while read id
do
 name=${id%%.*}
 echo "java -jar ~/biosoft/picard.jar MarkDuplicates I=${id} O=${name}.markdup.bam M=${name}.markdup.metrc.csv" >> markdup.commandlines
done

# ParaFly运行
ParaFly -c markdup.commandlines -CPU 16 2>markdup.err.log &
# rm -rf *.sorted.bam
# 上述这一步骤完成之后，需要对结果文件建立索引，以.bai结尾
for id in `ls *.markdup.bam`;do echo "java -jar ~/biosoft/picard.jar BuildBamIndex I=${id}" >> buildbamindex.commondlines;done
# ParaFly运行
ParaFly -c buildbamindex.commondlines -CPU 16 2>buildbamindex.err.log &

# 构建参考序列dict文件
java -jar ~/biosoft/picard.jar CreateSequenceDictionary R=../00.ref/arab_ref.fa O=../00.ref/arab_ref.dict
# 1、HaplotypeCaller
# mkdir 03.snpcalling && cd 03.snpcalling
ls ../02.bwa/*.bam | xargs -i basename {} > bam_id.txt
cat bam_id.txt | while read id
do
 name=${id%%.*}
 echo "gatk --java-options \"-Xmx15g\" HaplotypeCaller -R ../00.ref/arab_ref.fa -I ../02.bwa/$id -O ${name}.g.vcf.gz -ERC GVCF" >> haplotypecaller.commondlines
done

ParaFly -c haplotypecaller.commondlines -CPU 16 2>haplotypecaller.err.log &

# 2、GenomicsDBImport
ls *.g.vcf.gz > gvcf_id.txt
cut -d "." -f 1 gvcf_id.txt > sample_id.txt
paste -d "\t" sample_id.txt gvcf_id.txt > sample.map

my_database=arab_16_DB                    # 创建GenomicsDatabase的文件夹名称
samplemap=sample.map                      # 使用tab分隔符，第一列为sampleID，第二列为对应gvcf ID
interval_list=interval.list               # 每一行一个区域编号
gatk --java-options "-Xmx4g -Xms4g" \
GenomicsDBImport \
--genomicsdb-workspace-path $my_database \
--sample-name-map $samplemap \
-L $interval_list 2>genomicsDB_build.err.log &