/knowledge-of-xutengfei

群体遗传学

前言

我以往的群体遗传学工作主要是挖掘一个群体中基因型和表型之间的关联，外面很多公司都在做这类的研究，所以很多资料、文献都可以很容易地查到。

生物信息具体分析无非搞懂几件事情：

  1. 数据格式
  2. 基础统计学
  3. 软件参数意义
  4. 可视化
  5. 结果解读

数据格式

fasta 是最简单的序列格式。一条序列由两行组成，一行为以>开头的序列名称，一行为由ATCGN组成的序列。一般序列长度为100-200bp。

fastq 是测序仪下机序列的基本格式。一条序列由四行组成：序列名称、序列、+ 号、质量值。如果测序类型为双端测序（pair-end），则下机数据一个样有一对儿文件，通常文件以_1.fq.gz 和_2.fq.gz结尾。

gff3 是基因组注释文件的格式。其保存了一个基因在基因组上的具体区域，如果由 ncbi 或者 Ensembl 下载的注释文件，则还会包含基因具体的名称，细节自己查询。

vcf 是变体调用格式。几乎所有的变异类型（snp，indel，sv，cnv）都可以用 vcf 表示。vcf 由三部分组成：以##开头的变异元数据，包含了染色体信息，格式中字段的具体含义等。以#开头的变异内容头。变异具体内容，包含了变异位置、变异类型、具体每个样本的碱基类型等。

sam 序列比对生成文件的格式。包含了fastq文件的所有信息，还有其他比对信息，所以文件巨大，建议转成bam格式，项目结束后，删除比对文件。细节自己查询。

基础统计学

需要学习各个分布（二项，正态，泊松）的原理和应用，统计学检验的计算，特征相关性计算，箱线图运用。要懂得做研究的统计学前提：

样本数据够不够？（一组最少30个）

表型是不是正态分布？

缺失值是删除还是填补？

如果无法保证统计学前提，那么研究即使有显著的结果，那也是不可信的，片面的，无法应用到实践的。

软件参数意义

需要了解软件常用参数会对结果造成什么影响，一般来说，很基础的参数，软件作者会直接给最好的默认值，但是当问题出现的时候，第一个想到是不是自己操作有问题，第二个输入数据格式是不是有问题，最后就是软件参数是否有问题。

可视化

散点图：体现数据的分布

折线图：体现数据的变化趋势

条形图：不同组数据大小比较

饼图：各个成分之间比例

箱线图： 体现不同组数据直观分布，离群值检测

平行坐标图：不同组数据多维特征数值分布

结果解读

高情商：结果解读，低情商：讲故事

数据预处理

数据预处理就是去掉下机序列数据中不应该存在的碱基，或者质量不行的碱基。（接头、碱基N，低质量碱基） 我比较常用软件：cutadapt和TrimGalore，TrimGalore其实是cutadapt的套壳。如果不提供接头序列，他们都可以自己去统计，然后找一个最可能是接头的接头，然后再去除。

比对

序列比对通常用 bwa，现在bwa有了第二个版本，比对速度上有较大的提升，参数没有变化。https://github.com/bwa-mem2/bwa-mem2

snp indel

变异检测只用GATK，集群版本有v3.3 v3.7 v3.8 v4.x,记得对应java1.7和java1.8

变异过滤

变异过滤用vcftools，在几个方面对变异进行过滤：变异缺失率（90%），次等位基因频率maf（0.05），等位基因数（2)。网上可能对MAF解释有误导。

MAF解释: 比如一个变异位点为chr1-12345,变异类型为A-T,你做的研究共有100个样本（二倍体），那么这个位点就有200个碱基信息。碱基A叫做REF，T叫做ALT，T数量占200个碱基的比例就叫做次等位基因频率，只有当T数量占所有碱基的5%以上，这个位点才得以保留。当然现在有研究稀有变异的（MAF<0.05），但是不在今天讨论范畴。

cnv sv

拷贝数变异使用软件cnvcaller,结构变异使用软件breakdancer,smoove,manta。由于结构变异无法很好地确定准确变异边界，所以要用多个软件求结构并集。按以往实践来看manta最准确。

群体结构

群体结构分析使用软件admixture ，首先使用plink将变异文件转换为bed格式，然后使用admixture 进行计算。没有什么特别的参数，只有一个k参数注意，一般从2到10都跑一遍，意思是推测群体分2到10层。群体结构分层会影响后续的GWAS分析，因为GWAS的计算基础是线性模型，如果样本之间有关联，则违背了线性模型的统计学假设，所以要将群体G矩阵加入线性模型中（协变量）。

PCA

主成分分析和群体结构分析类似，也是寻找群体分层的方法，使用软件gcta可以很快的得出pca结果。参考：https://flystar233.github.io/2019/11/19/pca/

流程

https://github.com/flystar233/rstool

此流程是我从19年开始编写，维护2年的重测序流程，已经在多个项目中使用，如今版本号为：v1.6.5。如需要使用，还请发paper的时候引用 doi：10.5281/zenodo.4466539

GWAS

全基因组关联分析是将表型和基因型关联到一起的一种方法，基于线性模型。当表型是二分类时，使用逻辑回归模型，当表型是连续型时，使用一般线性模型或者混合线性模型。混合线性模型相比于一般线性模型可以添加 随机效应，随机效应在这里指代的是样本之间的亲缘关系（K矩阵）。使用软件EMMAX可以轻松计算变异位点显著性，其附属软件emmax-kin可以计算K矩阵。

显著位点

以以往经验，动植物物种变异信息经过严格的过滤后，进行gwas分析的位点大约在3002万个左右。gwas分析之后会给出每个位点的p值，p值越小，则和表型的关联显著性越高。那么如何定量到一个阈值呢？

  1. 最严格的Bonferroni 矫正。0.05/位点数 就是阈值。一般来说，这个数值都会达到10的负8次方，很少有位点的p值低到这个程度，从而无点可选。
  2. 计算有效位点数。软件simpleM_Ex通过输入变异文件，输入有效位点数，再使用0.05/有效位点数可以让阈值稍微降低。

变异注释

https://flystar233.github.io/2020/05/20/annovar/

基因组三代组装

基因组调查

使用软件jellyfish进行kmer分布计算，然后使用genomescope预估基因组大小，杂合性等基因组信息。

组装

三代组装软件繁多，只推荐几个好用的。

  1. shasta三代单装，几乎不用动参数
  2. wengan 二三代混装，准确性高
  3. raven 三代单装，contig长，自带两轮polish，由racon作者编写 组装成 contig 之后需要使用racon进行2轮三代纠错，使用pilon进行2-3轮的二代数据纠错，最后使用busco进行组装评估。如果不进行hic组装的话，可以用quickmerge软件综合2种组装结果，提高contig水平。

HIC组装

第一步

hic-pro 处理原始hic数据，得到有效片段。（没有改动）

第二步

python 重写getread.pl。此脚本是将 hic-pro结果有效片段文件中的所有序列名与原始 reads 匹配并提取出来。

与原代码相比，在同样的运行时间下，增加了对 fastq 文件新旧 2 种 ID 的兼容和程序进度提示。fastq 文件新 ID 示例:@V300075059L2C001R0020000012/1 ，旧 ID：@A00838:157:H2Y5TDSXY:4:1101:1398:1000 2:N:0:ATTCAGAA+CCTATCCT

使用方法
sh getread.sh rawdata_allValidPairs CL10_1.fq.gz CL10_2.fq.gz #建议qsub提交任务

第三步

3d-dna 本就是一套用 sh、perl、python、awk 编写的流程。

在最新版本中，3d-dna 所依赖的 juicer 开始使用 bwa 双端比对，3d-dna 不再要求输入预估的染色体数目。在原接口流程中（由 liuqun@genomics.cn 编写），程序一直会监测 juicer 生成的文件，如果程序开始生成 inter.hic，则 qdel 任务，并开始运行3d-dna。这样就存在了如果不是nohup提交任务，如果关掉 ssh 集群连接，则任务将会失败。

在原流程的基础上修改了代码，以适应新 juicer+3d-dna 参数，并直接生成 juicer 和 3d-dna程序所需代码， 当juicer 生成 merged_nodups.txt 文件后（即开始生成 inter.hic 时），可以杀死 juicer任务，开始运行 3d-dna。在运行过程中，可以把 0.hic(1.hic,2.hic,polished.hic 等等hic文件)和 0.asm(1.asm,2,asm 等等asm文件)导出到 juicerbox 中查看具体情况，通常情况下 0.hic 会挂载最多的原始基因组序列，rawchrom.hic 会确定最好的染色体边界。选择一个进行人工 review 热图，确定染色体边界之后导出 review.assembly 文件。然后再到集群运行 3d-dna review 程序得到最终的hic组装结果。

使用方法
python juicer_3d_dna_pipe.py --fq *.gz --ref refer.fa --res_name MboI --n_cpus 8
bash run-asm-pipeline-post-review.sh --sort-output  -r test.review.assembly contig.fa  ./merged_nodups.txt

脚本存放

/ldfssz1/MS_OP/USER/xutengfei1/script_py/hic_scr

/zfssz3/NASCT_/MS_PMO2017/xutengfei1/software/3d-dna_bin

具体参考 /ldfssz1/MS_OP/USER/xutengfei1/tangyu/new_HIC

编程

建议直接学习 python ，抛弃perl这种魔幻语言，即使你不干生信了，python也很有用处。 linux 中 awk sed sort grep cut 命令要会灵活应用。学习 python 基础就不用说了，主要在几个方面要精通

  1. 写代码清晰的方法（def）
  2. 学着写类（class）
  3. 写灵活的命令行接口（click）
  4. 文件、路径操作（pathlib）
  5. 脑子里必须记得每个正则符号代表什么，灵活组合（re）
  6. 建议把常用软件的读取，写入，格式转换，自己都重写一次，感觉真的不一样
  7. 各种常用的数据处理包的使用（csv，pandas，scikit-learn，pyod，jellyfish，pyfastx）

软件

linux

  1. fzf 模糊历史命令查询，你会回来感谢我的 https://github.com/junegunn/fzf
  2. ripgrep 类似 gerp，但是更强大 https://github.com/BurntSushi/ripgrep
  3. seqkit fasta、fastq 文件处理 https://github.com/shenwei356/seqkit
  4. csvtk 结构数据文件处理 https://github.com/shenwei356/csvtk
  5. miniconda3 编程、生信软件安装

windows

  1. sublime text3 文本编辑器
  2. cmder 更好的 cmd
  3. Termius ssh 连接集群（美观）
  4. Visual Studio Code 我一般用它方便debug
  5. moeditor 写 markdown文档用
  6. TED Notepad 比自带的文档编辑器好用
  7. pencil 流程图制作
  8. 多去 github 玩