生成的文件說明：

1.~/Process/fastqc_result 用于存放質(zhì)量評估的結(jié)果
2.~/Process/trim_out 用于存放過濾的結(jié)果
3.~/Process/result/depth.txt 將測序深度保存為depth.txt
4.~/Process/result/T_result.txt 將比對結(jié)果保存為T_result.txt
5.變異注釋生成avinput.variant_function和avinput.exonic_variant_function后綴的兩個結(jié)果文件保存在~/Process/result

1.安裝相關(guān)軟件及環(huán)境的腳本

#!/bin/bash
cd ~
mkdir ~/BioSofts
cd BioSofts
echo "..................................................................................................."
echo "配置java環(huán)境......"
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 "創(chuàng)建目錄用于安裝相關(guān)軟件......"
cd ~/BioSofts
echo ".................................................................................................."
echo "安裝Aspera軟件，用于下載大的測序數(shù)據(jù)文件......"
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 "下載格式轉(zhuǎn)換工具包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 "安裝測序數(shù)據(jù)質(zhì)控軟件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 "安裝測序數(shù)據(jù)的數(shù)據(jù)過濾軟件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 "安裝格式轉(zhuǎn)換軟件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.數(shù)據(jù)下載的腳本

#!/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 anonftp@ftp-private.ncbi.nlm.nih.gov:/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 anonftp@ftp-private.ncbi.nlm.nih.gov:/genomes/all/GCA/001/877/055/GCA_001877055.1_ASM187705v1/GCA_001877055.1_ASM187705v1_genomic.fna.gz ~/Seqs/
echo "下載完成！"
echo ".................................................................................."

echo ".................................................................................."
echo "下載參考基因組的gff文件帆竹，用于注釋的自定義數(shù)據(jù)庫......"
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 "..................數(shù)據(jù)準備完畢...................................................."

3.重測序數(shù)據(jù)分析腳本

#!/bin/bash
cd ~
mkdir Process
cp ~/Seqs/* ~/Process
echo "一、將序列與參考基因組的比對："
echo "1.數(shù)據(jù)文件的格式轉(zhuǎn)換"
cd Process
mkdir result
fastq-dump SRR022534.sra

echo "2.質(zhì)量評估"
mkdir fastqc_result #用于存放質(zhì)量評估的結(jié)果
fastqc SRR022534.fastq
mv *.zip *.html fastqc_result/

echo "3.測序數(shù)據(jù)的數(shù)據(jù)過濾"
mkdir trim_out #用于存放過濾的結(jié)果
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.測序數(shù)據(jù)比對到參考基因組得到sam文件"
bwa mem GCA_001877055.1_ASM187705v1_genomic.fna trim_out/SRR011534_out.fastq.gz >bwa_mem_SRR011534.sam

echo "6.sam文件轉(zhuǎn)換為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.測試參考基因組每個位點或一段區(qū)域的測序深度"
samtools depth bwa_mem_SRR011534.sorted.bam >result/depth.txt #將測序深度保存為depth.txt

echo "10.統(tǒng)計比對結(jié)果"
samtools flagstat bwa_mem_SRR011534.sorted.bam > result/T_result.txt #將比對結(jié)果保存為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.自定義注釋數(shù)據(jù)庫"
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后綴的兩個結(jié)果文件"
annotate_variation.pl --geneanno --dbtype refGene --buildver 7055-genome_new  bwa_mem_SRR011534.snps.avinput ~/BioSofts/annovar/humandb/
mv *_function result
echo "....................................第三部分完成.............................................."

echo ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>重測序數(shù)據(jù)分析完成>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"

4.start.sh

#!/bin/bash

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

echo "說明：本次運行需要以下軟件和環(huán)境的支持绢陌，并且需要將軟件加到環(huán)境目錄下挨下，請確認是否已安裝這些軟件和環(huán)境，如果是請輸入1,否請輸入0"
echo "Aspera sratoolkit FastQC Trimmomatic bwa samtools Annovar"
read response
if [ $response -eq 1  ] ; then
    echo "開始下載數(shù)據(jù)"
    exec ./data_download.sh
    echo "開始分析"
    exec ./precess.sh
elif [ $response -eq 0 ] ; then
    echo "開始安裝所需要的軟件"
    exec ./software.sh
    echo "開始下載數(shù)據(jù)"
    exec ./data_download.sh
    echo "開始分析"
    exec ./precess.sh
else
    echo "輸入有誤脐湾，退出臭笆！"
    exit 1
fi