已知達(dá)松維爾擬諾卡氏菌亞種是會導(dǎo)致人類放線菌瘤的環(huán)境生物郁稍，該樣本是來自Keddieii血根桿菌DSM 10542的通用樣品辜梳，旨通過基因組重測序探索其和參考基因組有何不同粱甫，找出基因組變異信息。

1.需要的軟件

? 軟件名：Aspera 版本號：4.0.2.38
? 軟件名：sratoolkit 版本號：2.11.1
? 軟件名：FastQC 版本號：0.11.7
? 軟件名：Trimmomatic版本號：0.38
? 軟件名：bwa 版本號：0.7.17-r1188
? 軟件名：samtools 版本號：1.7
? 軟件名：Annovar 版本：$Date: 2019-10-24 00:05:27 -0400 (Thu, 24 Oct 2019)
? 軟件名：bcftools 版本：1.7

2.數(shù)據(jù)下載

2.1下載sra數(shù)據(jù)

1.在NCBI官網(wǎng)的搜索框輸入SRR022534作瞄，可以看到該菌的目前研究以及進(jìn)展茶宵，在最下面Runs點擊SRR022534稽物，在跳轉(zhuǎn)界面點擊Data access股毫，會看到一個網(wǎng)址坛悉，復(fù)制網(wǎng)址奥额；
2.服務(wù)器建re_seq文件夾，wget下載亞種基因組文件椒惨。

wget https://sra-pub-runodp.s3.amazonaws.com/sra/SRR022534/SRR022534

或者prefetch下載

prefetch SRR022534

也可以選擇ascp下載拣度，不過如果基因組文件不是特別大不推薦這種方式俊性，比較麻煩搪桂，需要到EBI-ENA數(shù)據(jù)庫賦值ftp地址透敌。

2.2.下載參考基因組序列和基因組注釋文件

網(wǎng)址：ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/001/877/055/GCA_001877055.1_ASM187705v1/GCA_001877055.1_ASM187705v1_genomic.fna.gz
1.進(jìn)入NCBI官網(wǎng)，輸入GCA_001877055踢械，在界面右邊點擊FTP directory for GenBank assembly酗电，在下一個界面復(fù)制參考基因組.fna.gz文件和基因組注釋文件.gff.gz；
2.wget下載内列。

wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/001/877/055/GCA_001877055.1_ASM187705v1/GCA_001877055.1_ASM187705v1_genomic.fna.gz

wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/001/877/055/GCA_001877055.1_ASM187705v1/GCA_001877055.1_ASM187705v1_genomic.gff.gz

3.主要分析步驟和結(jié)果

3.1系列與參考基因組的比對

3.1.1數(shù)據(jù)文件的格式轉(zhuǎn)換

在NCBI數(shù)據(jù)庫檢索序列號SRR033534查看experiment得知測試數(shù)據(jù)類型為454的單末端測序撵术。

fastq-dump SRR022534.sra

3.1.2質(zhì)量評估

1.fastqc質(zhì)量評估，文件最好用絕對路徑话瞧；

touch fastp.sh
vim fastp.sh
#!/bin/bash
for i in 1 2 3 4
do
fastp -i SRR137492${i}.fastq.gz -o SRR137492${i}.fq.gz
done
#Esc : wq! 保存并退出
chmod 777 fastp.sh
./fastp.sh

2.本地查看fastqc report嫩与，由結(jié)果可知寝姿，測序數(shù)據(jù)質(zhì)量不是特別高，Per base sequence quality中有很多下四分位小于10或中位數(shù)小于25划滋，所以下一步要先數(shù)據(jù)過濾会油，去除低質(zhì)量序列。

3.1.3測序數(shù)據(jù)過濾

Trimmomatic參數(shù)說明

PE：雙端測序
-phread33：設(shè)置堿基的質(zhì)量格式
ILLUMINACLIP：
TruSeq3-PE.fa： <fastaWithAdaptersEtc>是fasta格式的接頭序列文件路徑
2： <seed mismatches>是將接頭序列的一段(不超過16bp)作為seed古毛，與reads進(jìn)行比對，能夠容忍的最大錯配(mismatch)數(shù)都许，這里是最多2個錯配
30： <palindrome clip threshold>是 a-R1, a-R2的比對分值閾值稻薇，達(dá)到閾值，才進(jìn)行切除胶征，這里設(shè)置閾值為30（大約比對50bp堿基）
10：<simple clip threshold>是任意(接頭)序列與read比對最低分閾值塞椎，大于這個閾值，才進(jìn)行切除睛低，這里設(shè)置閾值為10（大約比對17bp堿基）
LIDINGWINDOW:5:20：設(shè)置滑動窗口閾值案狠，以為5bp為窗口，這5bp堿基的平均質(zhì)量值低于20钱雷，要進(jìn)行切除
:LEADING:20：設(shè)置從reads起始開始骂铁，去除質(zhì)量低于閾值或為’N’的堿基，直到達(dá)到閾值不再去除罩抗，這里設(shè)置閾值為20
TRAILING:20：設(shè)置從reads末尾開始拉庵，去除質(zhì)量低于閾值的堿基或為’N’的堿基直到達(dá)到閾值不再去除，這里設(shè)置閾值為3套蒂，這種方法是去除特定的illumina平臺低質(zhì)量區(qū)域(由于illumina會將低質(zhì)量堿基標(biāo)記為2) 
MINLEN:75 設(shè)置read切除后的最短長度钞支，這里設(shè)置長度至少為36bp，長度小于36bp的reads被去除

fastqc /disk1/202031107010114/re_seq/trim_out/SRR022534.fastq.gz

3.1.4fastqc重新質(zhì)量評估

fastqc /disk1/202031107010114/re_seq/trim_out/SRR022534.fastq.gz

3.1.5建立參考基因組索引

gunzip disk1/202031107010114/re_seq/GCA_001877055.1_ASM187705v1_genomic.fna.gz
bwa index disk1/202031107010114/re_seq/GCA_001877055.1_ASM187705v1_genomic.fna

3.1.6測序數(shù)據(jù)比對到參考基因組得到sam文件

bwa mem disk1/202031107010114/re_seq/GCA_001877055.1_ASM187705v1_genomic.fna disk1/202031107010114/re_seq/trim_out/SRR022534_out.fastq.gz >bwa_mem_SRR022534.sam

3.1.7sam文件轉(zhuǎn)bam文件

samtools faidx disk1/202031107010114/re_seq/GCA_001877055.1_ASM187705v1_genomic.fna
samtools view -bhS -t GCA_001877055.1_ASM187705v1_genomic.fna.fai -o bwa_mem_SRR022534.bam bwa_mem_SRR022534.sam

3.1.8為bam文件排序

samtools sort bwa_mem_SRR022534.bam -o bwa_mem_SRR022534.sorted.bam

3.1.9為bam文件建立索引顯示基因組比對情況

samtools index bwa_mem_SRR022534.sorted.bam
samtools index bwa_mem_SRR022534.sorted.bam

3.1.10測試每個基因組每個位點的測試深度并統(tǒng)計比對結(jié)果

samtools depth bwa_mem_SRR022534.sorted.bam >>depth.txt
samtools flagstat bwa_mem_SRR022534.sorted.bam

3.2變異點檢測

3.2.1去除pcr重復(fù)

samtools rmdup bwa_mem_SRR022534.sorted.bam bwa_mem_SRR022534_nopcr.bam

3.2.2生成bcf文件

samtools mpileup -gf GCA_001877055.1_ASM187705v1_genomic.fna bwa_mem_SRR022534_nopcr.bam >bwa_mem_SRR022534.bcf

3.2.3基因變異檢測和篩選

bcftools call -vm bwa_mem_SRR022534.bcf -o bwa_mem_SRR022534.variants.bcf
bcftools view -v snps,indels bwa_mem_SRR022534.variants.bcf > bwa_mem_SRR022534.snps.vcf

3.3變異位點注釋

3.3.1生成annovar輸入文件

convert2annovar.pl -format vcf4 bwa_mem_SRR022534.snps.vcf > bwa_mem_SRR022534.snps.avinput

3.3.2自定義注釋數(shù)據(jù)框

gunzip disk1/202031107010114/re_seq/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_
cp 7055-genome_new_refGene* ~/Biosofts/annovar/humandb/

3.3.3注釋變異基因位點

annotate_variation.pl --geneanno --dbtype refGene --buildver 7055-genome_new  bwa_mem_SRR022534.snps.avinput ~/Biosofts/annovar/humandb/

4.批量處理的腳本

因為我所分析的sra文件只有一個操刀，無法更好地說明批量處理過程烁挟，就找了一組RNA-Seq數(shù)據(jù)，二者前面的分析流程基本一致骨坑。數(shù)據(jù)來自研究：抗體富集前后培養(yǎng)物中伯氏疏螺旋體基因表達(dá)中富含aBa的體外培養(yǎng)Bb代表SRR22149531和體外培養(yǎng)的Bb代表SRR22149536撼嗓。

4.1prefetch批量下載sra文件

touch prefetch.sh
vim prefetch.sh
#!/bin/bash
for i in 1 6
do
prefetch SRR2214953${i}
done
#Esc : wq! 保存并退出
chmod 777 prefetch.sh
./prefetch.sh

4.2fastq-dump批量從sra文件提取fastq文件

touch fastq_dump.sh
vim fastq_dump.sh
#!/bin/bash
for i in SRR*
do
echo ${i}#打印輸出文件名
fastq-dump --gzip ${i}/${i}.sra
#可根據(jù)數(shù)據(jù)需要選擇--spilt-files/--spilt-3/--spilt-spot
done
#Esc : wq! 保存并退出
chmod 777 fastq_dump.sh
./fastq_dump.sh

4.3fastqc批量質(zhì)控

touch fastqc.sh
vim fastqc.sh
#！/bin/bash
for i in 1 6
do
fastqc SRR2214953${i}
done
#Esc : wq! 保存并退出
chmod 777 fastqc.sh

4.4bwa批量比對并生成排好序的bam文件

touch bwa_samtools.sh
vim bwa_samtools.sh
#!/bin/bash
for i in 1 6;do
{
bwa index GCA_000181575.2_ASM18157v2_genomic.fna;
bwa mem GCA_000181575.2_ASM18157v2_genomic.fna SRR2214953${i}.fastq.gz >bwa_mem_SRR2214953${i}.sam;
samtools faidx GCA_000181575.2_ASM18157v2_genomic.fna;
samtools view -bhS -t GCA_000181575.2_ASM18157v2_genomic.fna.fai -o bwa_mem_SRR2214953${i}.bam bwa_mem_SRR2214953${i}.sam;
samtools sort bwa_mem_SRR2214953${i}.bam -o bwa_mem_SRR2214953${i}.sorted.bam;}
done
#Esc : wq! 保存并退出
chmod 777 bwa_samtools.sh
./bwa_samtools.sh

基因組重測序分析全過程

基因組重測序分析全過程

1.需要的軟件

2.數(shù)據(jù)下載

2.1下載sra數(shù)據(jù)

2.2.下載參考基因組序列和基因組注釋文件

3.主要分析步驟和結(jié)果

3.1系列與參考基因組的比對

3.1.1數(shù)據(jù)文件的格式轉(zhuǎn)換

3.1.2質(zhì)量評估

3.1.3測序數(shù)據(jù)過濾

3.1.4fastqc重新質(zhì)量評估

3.1.5建立參考基因組索引

3.1.6測序數(shù)據(jù)比對到參考基因組得到sam文件

3.1.7sam文件轉(zhuǎn)bam文件

3.1.8為bam文件排序

3.1.9為bam文件建立索引顯示基因組比對情況

3.1.10測試每個基因組每個位點的測試深度并統(tǒng)計比對結(jié)果

3.2變異點檢測

3.2.1去除pcr重復(fù)

3.2.2生成bcf文件

3.2.3基因變異檢測和篩選

3.3變異位點注釋

3.3.1生成annovar輸入文件

3.3.2自定義注釋數(shù)據(jù)框

3.3.3注釋變異基因位點

4.批量處理的腳本

4.1prefetch批量下載sra文件

4.2fastq-dump批量從sra文件提取fastq文件

4.3fastqc批量質(zhì)控

4.4bwa批量比對并生成排好序的bam文件

基因組重測序分析全過程

1.需要的軟件

2.數(shù)據(jù)下載

2.1下載sra數(shù)據(jù)

2.2.下載參考基因組序列和基因組注釋文件

3.主要分析步驟和結(jié)果

3.1系列與參考基因組的比對

3.1.1數(shù)據(jù)文件的格式轉(zhuǎn)換

3.1.2質(zhì)量評估

3.1.3測序數(shù)據(jù)過濾

3.1.4fastqc重新質(zhì)量評估

3.1.5建立參考基因組索引

3.1.6測序數(shù)據(jù)比對到參考基因組得到sam文件

3.1.7sam文件轉(zhuǎn)bam文件

3.1.8為bam文件排序

3.1.9為bam文件建立索引 顯示基因組比對情況

3.1.10測試每個基因組每個位點的測試深度并統(tǒng)計比對結(jié)果

3.2變異點檢測

3.2.1去除pcr重復(fù)

3.2.2生成bcf文件

3.2.3基因變異檢測和篩選

3.3變異位點注釋

3.3.1生成annovar輸入文件

3.3.2自定義注釋數(shù)據(jù)框

3.3.3注釋變異基因位點

4.批量處理的腳本

4.1prefetch批量下載sra文件

4.2fastq-dump批量從sra文件提取fastq文件

4.3fastqc批量質(zhì)控

4.4bwa批量比對并生成排好序的bam文件

3.1.9為bam文件建立索引顯示基因組比對情況