基因组测序数据分析linux,基因组重测序数据分析脚本
生成的文件说明：
1.~/Process/fastqc_result 用于存放质量评估的结果
2.~/Process/trim_out 用于存放过滤的结果
3.~/Process/result/depth.txt 将测序深度保存为depth.txt
4.~/Process/result/T_result.txt 将比对结果保存为T_result.txt
在~/Process/result目录下创建并保存avinput.variant_function和avinput.exonic_variant_function两个结果文件
1.安装相关软件及环境的脚本
#!/bin/bash
cd ~
mkdir ~/BioSofts
cd BioSofts
echo "..................................................................................................."
echo "配置java环境......"
sudo mkdir /usr/java
使用wget命令下载JDK 8u131版本的Linux-x64系统包，请注意以下参数设置：禁止证书验证（skip certificate validation）、不读取Cookies（without cookies），并附加HTTP头信息（HTTP header）设置Cookie为accept-securebackup-cookie。
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函数将CLASSPATH变量的值以特定路径拼接的方式赋值给文件/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......"
bash从美国达利实验室的网站上下载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......"
Use wget to obtain the file from the specified position: 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 "..................................................................................."
获取达松维尔拟诺卡氏菌亚种（cardiopsis dassonvillei subsp.）的相关基因组数据
ascp --transpose --input-file .aspera(asepa ID DSA SSH config file path等价于openssh) --key-limit 1 --max-file-length 200M anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR/SRR022/SRR/SRR022534/SRR022534.sra ~/personal directory path
echo "下载完成!"
echo "..................................................................................."
echo "..................................................................................."
echo "达松维尔拟诺卡氏菌(Nocardiopsis dassonvillei )参考基因组序列......"
使用ascp工具执行同步操作，并连接到指定的SSH地址,设置同步选项为端口数为1,本地端口范围从20到2万,目标用户的IP地址为anonftp@ftp-private.ncbi.nlm.nih.gov:/genomes/all/GCA/...,并传输指定的基因组数据至本地路径~/Seqs.
echo "下载完成！"
echo ".................................................................................."
echo ".................................................................................."
echo "下载参考基因组的gff文件，用于注释的自定义数据库......"
通过使用wget命令从FTP服务器访问并下载指定的基因组数据文件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命令行工具运行TrimMatic版本号为...的JAR文件，在特定目录下执行 trimming操作。
该脚本通过Phred score 99进行质量控制过滤，并对输入的fastq序列进行处理。
将ILLUMINACLIP配置为...序列 adapters文件夹下的...序列，并根据所需参数设置剪切策略。
具体来说：
ILLUMINACLIP配置为...序列 adapters文件夹下的...序列；
将滑动窗口大小设为5；
设置leading和trailing剪切长度均为20；
设定最小长度阈值为75。
输出结果存储在指定的目标fastq.gz文件中。
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 examine -bhS -t GCA_001877055.1_ASM187705v1_genomic.fna.fai -output=bwa_mem_SRR011534.sam 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.统计比对结果"
该脚本会将比对结果输出到T_result.txt文件中
echo "......................................第一部分完成..............................................."
echo "二、变异位点的检测"
echo "1.去除PCR重复"
samtools rmdup bwa_mem_SRR011534.sorted.bam bwa_mem_SRR011534_nopcr.bam
echo "2.生成bcf文件"
SAMTOOLS PILEUP -GF GCA_001877055.2 Genomic.fna Bamboos Mem SRR XXXXXXXXXX Bamboos Mem SRR XXXXXXXXXX.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.变异位点的过滤"
Apply the bcftools filter with options --output bwa_mem_SRR011534.snps.filtered.vcf --from 'QUAL≥20 AND DP≥5' to process the input file 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
sh scripts/retrieve_seq_from_fasta.sh --format refGene using GCA_001877055.1_ASM187705v1_genomic.fna; output files: 7055-genome_new_refGeneMrna.fa and 7055-genome_new_refGene.txt
cp 7055-genome_new_refGene* ~/BioSofts/annovar/humandb/
标记变异基因位点，并创建具有 avinput.variant_function 和 avinput.exonic_variant_function 后缀的两个输出文件
script annotate_variation [--annot_type] [--reference_genome] [--reference_version 7055新版本] bwa_mem_SRR011634.snps.input ~/BioSofts/annovar/humandb/
mv *_function result
echo "....................................第三部分完成.............................................."
echo ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>重测序数据分析完成>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"
#!/bin/bash
chmod 777 software.sh
chmod 777 process.sh
chmod 777 data_download.sh
在运行时，请核实您是否已经安装了必要的软件和服务？如果您已经安装了这些软件和服务，请输入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