基因组重测序数据分析脚本

生成的文件说明：

1.~/Process/fastqc_result 用于存放质量评估的结果
2.~/Process/trim_out 用于存放过滤的结果
3.~/Process/result/depth.txt 将测序深度保存为depth.txt
4.~/Process/result/T_result.txt 将比对结果保存为T_result.txt
5.变异注释生成avinput.variant_function和avinput.exonic_variant_function后缀的两个结果文件保存在~/Process/result

1.安装相关软件及环境的脚本

#!/bin/bash
cd ~
mkdir ~/BioSofts
cd BioSofts
echo "..................................................................................................."
echo "配置java环境......"
sudo mkdir /usr/java
wget --no-check-certificate --no-cookies --header "Cookie: oraclelicense=accept-securebackup-cookie" http://download.oracle.com/otn-pub/java/jdk/8u131-b11/d54c1d3a095b4ff2b6607d096fa80163/jdk-8u131-linux-x64.tar.gz
sudo tar -zvxf jdk-8u131-linux-x64.tar.gz -C /usr/java/
cd /usr/java
sudo ln -s jdk1.8.0_131 latest
sudo ln -s /usr/java/latest default
sudo chmod 777 /etc/profile
sudo echo 'export JAVA_HOME=/usr/java/latest' >> /etc/profile
sudo echo 'export PATH=$JAVA_HOME/bin:$JAVA_HOME/jre/bin:$PATH' >> /etc/profile
sudo echo 'export CLASSPATH=/usr/java:$JAVA_HOME/lib/dt.jar:$JAVA_HOME/lib/tools.jar' >> /etc/profile
source /etc/profile
sudo chmod 440 /etc/profile
echo "配置完成！"
echo "..................................................................................................."

echo "..................................................................................................."
echo "创建目录用于安装相关软件......"
cd ~/BioSofts
echo ".................................................................................................."
echo "安装Aspera软件，用于下载大的测序数据文件......"
wget https://download.asperasoft.com/download/sw/connect/3.6.2/aspera-connect-3.6.2.117442-linux-64.tar.gz
tar zvxf aspera-connect-3.6.2.117442-linux-64.tar.gz
bash aspera-connect-3.6.2.117442-linux-64.sh
echo 'export PATH=$PATH:~/.aspera/connect/bin'>>~/.bashrc
source ~/.bashrc
echo "Aspera安装完成!"
echo "..................................................................................................."

echo "..................................................................................................."
echo "下载格式转换工具包sratoolkit......"
wget https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.9.2/sratoolkit.2.9.2-ubuntu64.tar.gz
tar -zvxf sratoolkit.2.9.2-ubuntu64.tar.gz -C ~/BioSofts/
echo 'export PATH=~/BioSofts/sratoolkit.2.9.2-ubuntu64/bin:$PATH' >> ~/.bashrc
source ~/.bashrc
echo "sratoolkit安装完成！"
echo "..................................................................................................."


echo "..................................................................................................."
echo "安装测序数据质控软件FastQc......"
cd ~/BioSofts
wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.7.zip
unzip fastqc_v0.11.7.zip
echo 'export PATH=~/BioSofts/FastQC:$PATH'>>~/.bashrc
source ~/.bashrc
chmod 777 ~/BioSofts/FastQC/fastqc
echo "安装完成！"
echo "..................................................................................................."

echo "..................................................................................................."
echo "安装测序数据的数据过滤软件Trimmomatic......"
wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.38.zip
unzip Trimmomatic-0.38.zip
echo "安装完成！"
echo "..................................................................................................."

echo "..................................................................................................."
echo "安装比对软件bwa......"
sudo apt install bwa
echo "安装完成！"
echo "..................................................................................................."

echo "..................................................................................................."
echo "安装格式转换软件samtools......"
sudo apt install samtools
echo "安装完成！"
echo "..................................................................................................."

echo "..................................................................................................."
echo "变异基因注释软件Annovar......"
wget http://www.openbioinformatics.org/annovar/download/0wgxR2rIVP/annovar.latest.tar.gz
tar -zvxf annovar.latest.tar.gz
echo 'export PATH=~/BioSofts/annovar:$PATH'>>~/.bashrc
source ~/.bashrc
echo "..................................................................................................."

2.数据下载的脚本

#!/bin/bash 
cd ~
mkdir Seqs
cd Seqs
echo "..................................................................................."
echo "下载达松维尔拟诺卡氏菌亚种（cardiopsis dassonvillei subsp. dassonvillei）的基因组......"
ascp -T -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh -k 1 -l 200m [email protected]:/sra/sra-instant/reads/ByRun/sra/SRR/SRR022/SRR022534/SRR022534.sra ~/Seqs/
echo "下载完成!"
echo "..................................................................................."

echo "..................................................................................."
echo "达松维尔拟诺卡氏菌（Nocardiopsis dassonvillei ）参考基因组序列......"
ascp -T -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh -k 1 -l 200m [email protected]:/genomes/all/GCA/001/877/055/GCA_001877055.1_ASM187705v1/GCA_001877055.1_ASM187705v1_genomic.fna.gz ~/Seqs/
echo "下载完成！"
echo ".................................................................................."

echo ".................................................................................."
echo "下载参考基因组的gff文件，用于注释的自定义数据库......"
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/001/877/055/GCA_001877055.1_ASM187705v1/GCA_001877055.1_ASM187705v1_genomic.gff.gz
echo "下载完成！"
echo "..................数据准备完毕...................................................."

3.重测序数据分析脚本

#!/bin/bash
cd ~
mkdir Process
cp ~/Seqs/* ~/Process
echo "一、将序列与参考基因组的比对："
echo "1.数据文件的格式转换"
cd Process
mkdir result
fastq-dump SRR022534.sra

echo "2.质量评估"
mkdir fastqc_result #用于存放质量评估的结果
fastqc SRR022534.fastq
mv *.zip *.html fastqc_result/

echo "3.测序数据的数据过滤"
mkdir trim_out #用于存放过滤的结果
java -jar ~/BioSofts/Trimmomatic-0.38/trimmomatic-0.38.jar SE -phred33 SRR022534.fastq ./trim_out/SRR011534_out.fastq.gz ILLUMINACLIP:/home/lizeguo/BioSofts/Trimmomatic-0.38/adapters/TruSeq2-PE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:20 TRAILING:20 MINLEN:75

echo "4.建立参考基因组索引"
gunzip GCA_001877055.1_ASM187705v1_genomic.fna.gz
bwa index GCA_001877055.1_ASM187705v1_genomic.fna

echo "5.测序数据比对到参考基因组得到sam文件"
bwa mem GCA_001877055.1_ASM187705v1_genomic.fna trim_out/SRR011534_out.fastq.gz >bwa_mem_SRR011534.sam

echo "6.sam文件转换为bam文件"
samtools faidx GCA_001877055.1_ASM187705v1_genomic.fna
samtools view -bhS -t GCA_001877055.1_ASM187705v1_genomic.fna.fai -o bwa_mem_SRR011534.bam bwa_mem_SRR011534.sam 

echo "7.为bam文件排序"
samtools sort bwa_mem_SRR011534.bam -o bwa_mem_SRR011534.sorted.bam

echo "8.为bam文件建立索引"
samtools index bwa_mem_SRR011534.sorted.bam

echo "9.测试参考基因组每个位点或一段区域的测序深度"
samtools depth bwa_mem_SRR011534.sorted.bam >result/depth.txt #将测序深度保存为depth.txt

echo "10.统计比对结果"
samtools flagstat bwa_mem_SRR011534.sorted.bam > result/T_result.txt #将比对结果保存为T_result.txt
echo "......................................第一部分完成..............................................."

echo "二、变异位点的检测"
echo "1.去除PCR重复"
samtools rmdup bwa_mem_SRR011534.sorted.bam bwa_mem_SRR011534_nopcr.bam

echo "2.生成bcf文件"
samtools mpileup -gf GCA_001877055.1_ASM187705v1_genomic.fna bwa_mem_SRR011534_nopcr.bam >bwa_mem_SRR011534.bcf

echo "3.基因变异检测，得到bwa_mem_SRR011534.snps.bcf文件"
bcftools call -vm bwa_mem_SRR011534.bcf -o bwa_mem_SRR011534.variants.bcf
bcftools view -v snps,indels bwa_mem_SRR011534.variants.bcf > bwa_mem_SRR011534.snps.vcf

echo "4.变异位点的过滤"
bcftools filter -o bwa_mem_SRR011534.snps.filtered.vcf -i 'QUAL>20 &&DP>5' bwa_mem_SRR011534.snps.vcf
echo ".....................................第二部分完成..............................................."

echo "三、变异基因注释"
echo "1.生成annovar输入文件"
convert2annovar.pl -format vcf4 bwa_mem_SRR011534.snps.vcf > bwa_mem_SRR011534.snps.avinput

echo "2.自定义注释数据库"
gunzip GCA_001877055.1_ASM187705v1_genomic.gff.gz
wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/gff3ToGenePred
chmod 777 gff3ToGenePred 
./gff3ToGenePred -useName GCA_001877055.1_ASM187705v1_genomic.gff 7055-genome_refGene.txt
cut -f 12 7055-genome_refGene.txt >column1.txt
cut -f 2-15 7055-genome_refGene.txt >column_else.txt
paste column1.txt column_else.txt >7055-genome_new_refGene.txt
retrieve_seq_from_fasta.pl -format refGene -seqfile GCA_001877055.1_ASM187705v1_genomic.fna -outfile 7055-genome_new_refGeneMrna.fa 7055-genome_new_refGene.txt
cp 7055-genome_new_refGene* ~/BioSofts/annovar/humandb/

echo "3.注释变异基因位点，生成avinput.variant_function和avinput.exonic_variant_function后缀的两个结果文件"
annotate_variation.pl --geneanno --dbtype refGene --buildver 7055-genome_new  bwa_mem_SRR011534.snps.avinput ~/BioSofts/annovar/humandb/
mv *_function result
echo "....................................第三部分完成.............................................."

echo ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>重测序数据分析完成>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"

4.start.sh

#!/bin/bash

chmod 777 software.sh
chmod 777 process.sh
chmod 777 data_download.sh

echo "说明：本次运行需要以下软件和环境的支持，并且需要将软件加到环境目录下，请确认是否已安装这些软件和环境，如果是请输入1,否请输入0"
echo "Aspera sratoolkit FastQC Trimmomatic bwa samtools Annovar"
read response
if [ $response -eq 1  ] ; then
    echo "开始下载数据"
    exec ./data_download.sh
    echo "开始分析"
    exec ./precess.sh
elif [ $response -eq 0 ] ; then
    echo "开始安装所需要的软件"
    exec ./software.sh
    echo "开始下载数据"
    exec ./data_download.sh
    echo "开始分析"
    exec ./precess.sh
else
    echo "输入有误，退出！"
    exit 1
fi