此博客是对列日大学Marc HANIKENNE的Practical Genomics 课程的总结， 并且我将其分为四个部分：
第一部分：测序数据的质控和过滤
第二部分：遗传variants的鉴定和表征
第三部分：将一个插入突变映射到基因组
第四部分：使用shell将上述三个步骤串联起来，变成pipeline，流水作业。
0 目标与数据
目标
- 使用bwa-mem, samtools, IGV, bcftools 和R代码将reads映射到参考基因组，并且鉴定出SNP和indel variants.
- 使用velvet和blast从头组装一株绿藻莱茵衣藻基因组，并且鉴定出其中的插入突变。
数据
使用真实数据集练习
三种衣藻菌株（WT1、WT2 和突变体）使用 NextSeq Illumina 测序仪以高通量模式（400 Mio. 簇）进行测序，以获得 75 bp 的配对末端读数。 在这里，我们将首先访问 Durandal 上的序列文件并确定每个菌株的读取计数。 衣藻基因组大小约为 120 Mb，原始测序深度是多少？
1 质控使用的是fastqc
fastqc *.fastq
查看html文件，对得到数据初步查看，得到过滤使用的参数
过滤使用：trimmomatic， 如下：
再次使用fastqc查看数据质量：
2 鉴定遗传variants
需要与参考基因组比对。练习中对SNPs和indels进行鉴定
read mapping ,将使用BWAS软件， 并且使用BWA-MEM算法，其能更快和准确的检索。详细查看：
http://bio-bwa.sourceforge.net/bwa.shtml
BWA-MEM结果输入为SAM格式的文件，方便人阅读。还有一种BAM格式的文件。是便于电脑阅读的二进制文件。
SAM的详细介绍：
https://samtools.github.io/hts-specs/
对于SAM和BAM数据的处理，使用samtools(http://www.htslib.org/doc/samtools.html)。
samtools和bcftools联合使用可以对BAM文件进行call和filter variancts。
bcftools（http://www.htslib.org/doc/bcftools.html）能够从VCF和BCF格式文件中调用和处理variants。
BCF/VCF格式可以查看：https://samtools.github.io/hts-specs/
variant调用是比较麻烦的问题，可以查看：
https://github.com/samtools/bcftools/wiki/HOWTOs#mpileup-calling
http://samtools.github.io/bcftools/howtos/index.html
在鉴定variants时， 您确实必须确保可以将真实的variants与错误和偏差区分开来。
以下将使用基础的bcftools软件进行分析。
2.1 下载参考基因组和注释文件
找对自己研究物种的参考基因组，下载fasta和对应的gff文件
2.2 对下载的genome建立引所（index）。
将使用bwas index:
使用代码：bwa index -p ref *fasta
会得到五个文件：
2.3 读取参考基因组的mapping
使用BWA-MEM算法，并且以samtools应对SAM
mapping:
SAM 转为BAM
对SAM排序，并且对其进行index
2.4 对于mapping data进行可视化
使用IGV（integrative genomics viewer）, download网址:
http://software.broadinstitute.org/software/igv/
用到的文件有：参考基因组的index， *bam, *.bam.bai文件（后两个都是已经排序）
2.5 variant calling
如果只查看和操作variant可能比查看全基因组要好一些，所以需要对BAM文件进行剔除和过滤其他的信息，只剩下variants,这个过程称为call variants.
使用samtools软件
为一个或多个 BAM 文件生成 VCF、BCF 或堆积。使用bcftools文件，
下载VCF格式的文件，可以在IGV软件进行查看veriants.
2.6 提取进行统计
使用bcftools的stats命令对vcf.gz进行提取
2.7 对统计结果可视化
使用
plot-vcfstats -p plots_wt1_mutant wt1_mutant_comp.stats
查看*-summary.pdf文件
2.8 画韦恩图
下载vcf文件，在解压
使用R画图，VennDiagram包
3 在基因组找到基因组突变
突变株已使用插入诱变获得：WT 株已被转化为带有潮霉素 (Hygr) 抗性基因的基因盒。 该盒随机插入基因组中。 基于 Hygr 抗性选择数千个转化体并筛选感兴趣的表型。 在这一部分中，我们将使用一种简单的方法来识别基因组中盒式磁带的插入位点，这是识别致病基因的第一步。
3.1 基因组组装
使用velvet进行
统计contigs数
grep \> assembly_mutant/contigs.fa
检查quast结果
less quast_results/latest/report.txt
cp assembly_mutant/contigs.fa mutant.fasta
3.2 scaffolding
使用SSPACE来提高scaffolding
查看scaffolds数量
grep -c > sspace.final.scaffolds.fasta
3.3 查看组装的状态
使用QUAST进行
3.4 建一个blast库
使用makeblastdb，其默认在Biolinux已下载。
https://www.ncbi.nlm.nih.gov/books/NBK279688/
3.5 Blast分析
命令行说明：https://www.ncbi.nlm.nih.gov/books/NBK279690/
blastn -help # 进行查看options
命令：
blastn -db <your_name>_db -query ../../hygr.txt -out <your_name>_bn -evalue 1 -num_alignments 1000 -outfmt 6
结果在 <your_name>_dn
格式和结果如下：
3.6 复制鉴定出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网站上进行比对
结果序列在：<your_name>_hygr_scaffolds.fasta， 将其上传到Phytozome, 选择对应的基因组，即可比对得出结果（找到突变的基因）。
