此博客是對列日大學(xué)Marc HANIKENNE的Practical Genomics 課程的總結(jié)， 并且我將其分為四個部分：
第一部分：測序數(shù)據(jù)的質(zhì)控和過濾
第二部分：遺傳variants的鑒定和表征
第三部分：將一個插入突變映射到基因組
第四部分：使用shell將上述三個步驟串聯(lián)起來玖媚，變成pipeline挠唆，流水作業(yè)。

0 目標(biāo)與數(shù)據(jù)

目標(biāo)

1. 使用bwa-mem, samtools, IGV, bcftools 和R代碼將reads映射到參考基因組雹仿，并且鑒定出SNP和indel variants.
2. 使用velvet和blast從頭組裝一株綠藻萊茵衣藻基因組增热，并且鑒定出其中的插入突變。

數(shù)據(jù)

使用真實(shí)數(shù)據(jù)集練習(xí)
三種衣藻菌株（WT1胧辽、WT2 和突變體）使用 NextSeq Illumina 測序儀以高通量模式（400 Mio. 簇）進(jìn)行測序峻仇，以獲得 75 bp 的配對末端讀數(shù)。 在這里邑商，我們將首先訪問 Durandal 上的序列文件并確定每個菌株的讀取計(jì)數(shù)摄咆。 衣藻基因組大小約為 120 Mb，原始測序深度是多少人断？

1 質(zhì)控使用的是fastqc

fastqc *.fastq

查看html文件吭从，對得到數(shù)據(jù)初步查看，得到過濾使用的參數(shù)
過濾使用：trimmomatic恶迈， 如下：

image.png

再次使用fastqc查看數(shù)據(jù)質(zhì)量：

2 鑒定遺傳variants

需要與參考基因組比對涩金。練習(xí)中對SNPs和indels進(jìn)行鑒定
read mapping ,將使用BWAS軟件， 并且使用BWA-MEM算法，其能更快和準(zhǔn)確的檢索步做。詳細(xì)查看：
http://bio-bwa.sourceforge.net/bwa.shtml

BWA-MEM結(jié)果輸入為SAM格式的文件副渴，方便人閱讀。還有一種BAM格式的文件全度。是便于電腦閱讀的二進(jìn)制文件煮剧。
SAM的詳細(xì)介紹：
https://samtools.github.io/hts-specs/

對于SAM和BAM數(shù)據(jù)的處理，使用samtools(http://www.htslib.org/doc/samtools.html)将鸵。
samtools和bcftools聯(lián)合使用可以對BAM文件進(jìn)行call和filter variancts勉盅。
bcftools（http://www.htslib.org/doc/bcftools.html）能夠從VCF和BCF格式文件中調(diào)用和處理variants。
BCF/VCF格式可以查看：https://samtools.github.io/hts-specs/

variant調(diào)用是比較麻煩的問題顶掉，可以查看：
https://github.com/samtools/bcftools/wiki/HOWTOs#mpileup-calling
http://samtools.github.io/bcftools/howtos/index.html

在鑒定variants時菇篡， 您確實(shí)必須確保可以將真實(shí)的variants與錯誤和偏差區(qū)分開來一喘。
以下將使用基礎(chǔ)的bcftools軟件進(jìn)行分析驱还。

2.1 下載參考基因組和注釋文件

找對自己研究物種的參考基因組，下載fasta和對應(yīng)的gff文件

2.2 對下載的genome建立引所（index）凸克。

將使用bwas index:

image.png

image.png

使用代碼：bwa index -p ref *fasta
會得到五個文件：

image.png

2.3 讀取參考基因組的mapping

使用BWA-MEM算法议蟆，并且以samtools應(yīng)對SAM
mapping:

image.png

SAM 轉(zhuǎn)為BAM

image.png

對SAM排序，并且對其進(jìn)行index

image.png

2.4 對于mapping data進(jìn)行可視化

使用IGV（integrative genomics viewer）, download網(wǎng)址:
http://software.broadinstitute.org/software/igv/
用到的文件有：參考基因組的index萎战， *bam, *.bam.bai文件（后兩個都是已經(jīng)排序）

2.5 variant calling

如果只查看和操作variant可能比查看全基因組要好一些咐容，所以需要對BAM文件進(jìn)行剔除和過濾其他的信息，只剩下variants,這個過程稱為call variants.
使用samtools軟件

image.png

為一個或多個 BAM 文件生成 VCF蚂维、BCF 或堆積戳粒。使用bcftools文件，

image.png

下載VCF格式的文件虫啥，可以在IGV軟件進(jìn)行查看veriants.

2.6 提取進(jìn)行統(tǒng)計(jì)

使用bcftools的stats命令對vcf.gz進(jìn)行提取

image.png

2.7 對統(tǒng)計(jì)結(jié)果可視化

使用

plot-vcfstats -p plots_wt1_mutant wt1_mutant_comp.stats

查看*-summary.pdf文件

2.8 畫韋恩圖

下載vcf文件蔚约，在解壓
使用R畫圖，VennDiagram包

3 在基因組找到基因組突變

突變株已使用插入誘變獲得：WT 株已被轉(zhuǎn)化為帶有潮霉素 (Hygr) 抗性基因的基因盒涂籽。 該盒隨機(jī)插入基因組中苹祟。 基于 Hygr 抗性選擇數(shù)千個轉(zhuǎn)化體并篩選感興趣的表型。 在這一部分中评雌，我們將使用一種簡單的方法來識別基因組中盒式磁帶的插入位點(diǎn)树枫，這是識別致病基因的第一步。

3.1 基因組組裝

使用velvet進(jìn)行

統(tǒng)計(jì)contigs數(shù)

grep \> assembly_mutant/contigs.fa

檢查quast結(jié)果

less quast_results/latest/report.txt

cp assembly_mutant/contigs.fa mutant.fasta

3.2 scaffolding

使用SSPACE來提高scaffolding

查看scaffolds數(shù)量
grep -c > sspace.final.scaffolds.fasta

3.3 查看組裝的狀態(tài)

使用QUAST進(jìn)行

3.4 建一個blast庫

使用makeblastdb景东，其默認(rèn)在Biolinux已下載砂轻。
https://www.ncbi.nlm.nih.gov/books/NBK279688/

image.png

3.5 Blast分析

命令行說明：https://www.ncbi.nlm.nih.gov/books/NBK279690/
blastn -help # 進(jìn)行查看options
命令：

blastn -db <your_name>_db -query ../../hygr.txt -out <your_name>_bn -evalue 1 -num_alignments 1000 -outfmt 6

結(jié)果在 <your_name>_dn
格式和結(jié)果如下：

image.png

image.png

3.6 復(fù)制鑒定出scaffolds的序列

使用blast工具中的blastdbcmd

代碼：

cut -f2 <your_name>_bn > <your_name>_scaffolds.ids
blastdbcmd -db <your_name>_db -entry_batch <your_name>_scaffolds.ids >  <your_name>_hygr_scaffolds.fasta

3.7 在Phytozome網(wǎng)站上進(jìn)行比對

結(jié)果序列在：<your_name>_hygr_scaffolds.fasta， 將其上傳到Phytozome, 選擇對應(yīng)的基因組斤吐，即可比對得出結(jié)果（找到突變的基因）搔涝。

4 實(shí)際分析中丝里，可以寫一個shell命令文件，請上述步驟一次完成体谒，實(shí)現(xiàn)分析的自動化和批量進(jìn)行