❶ ATAC-seq专题---生信分析流程
ATAC-seq信息分析流程主要分为以下几个部分：数据质控、序列比对、峰检测、motif分析、峰注释、富集分析，下面将对各部分内容进行展开讲解。
下机数据经过过滤去除接头含量过高或低质量的reads，得到clean reads用于后续分析。常见的trim软件有Trimmomatic、Skewer、fastp等。fastp是一款比较新的软件，使用时可以用--adapter_sequence/--adapter_sequence_r2参数传入接头序列，也可以不填这两个参数，软件会自动识别接头并进行剪切。如：
fastp \
--in1 A1_1.fq.gz \ # read1原始fq文件
--out1 A1_clean_1.fq.gz \ # read1过滤后输出的fq文件
--in2 A1_2.fq.gz \ # read2原始fq文件
--out2 A1_clean_2.fq.gz \ # read2过滤后输出的fq文件
--cut_tail \ #从3’端向5’端滑窗，如果窗口内碱基的平均质量值小于设定阈值，则剪切
--cut_tail_window_size=1 \ #窗口大小
--cut_tail_mean_quality=30 \ #cut_tail参数对应的平均质量阈值
--average_qual=30 \ #如果一条read的碱基平均质量值小于该值即会被舍弃
--length_required=20 \ #经过剪切后的reads长度如果小于该值会被舍弃
fastp软件的详细使用方法可参考：https://github.com/OpenGene/fastp。fastp软件对于trim结果会生成网页版的报告，可参考官网示例http://opengene.org/fastp/fastp.html和http://opengene.org/fastp/fastp.json，也可以用FastQC软件对trim前后的数据质量进行评估，FastQC软件会对单端的数据给出结果，如果是PE测序需要分别运行两次来评估read1和read2的数据质量。
如：
fastqc A1_1.fq.gz
fastqc A1_2.fq.gz
FastQC会对reads从碱基质量、接头含量、N含量、高重复序列等多个方面对reads质量进行评估，生成详细的网页版报告，可参考官网示例：http://www.bioinformatics.babraham.ac.uk/projects/fastqc/good_sequence_short_fastqc.html
经过trim得到的reads可以使用BWA、bowtie2等软件进行比对。首先需要确定参考基因组fa文件，对fa文件建立索引。不同的软件有各自建立索引的命令，BWA软件可以参考如下方式建立索引：
bwa index genome.fa
建立好索引后即可开始比对，ATAC-seq推荐使用mem算法，输出文件经samtools排序输出bam：
bwa mem genome.fa A1_clean_1.fq.gz A1_clean_2.fq.gz
| samtools sort -O bam -T A1 > A1.bam
值得注意的是，在实验过程中质体并不能完全去除，因此会有部分reads比对到质体序列上，需要去除比对到质体上的序列，去除质体序列可以通过samtools提取，具体方法如下：首先将不含质体的染色体名称写到一个chrlist文件中，一条染色体的名称写成一行，然后执行如下命令即可得到去除质体的bam
samtools view -b A1.bam $chrlist > A1.del_MT_PT.bam
用于后续分析的reads需要时唯一比对且去重复的，bwa比对结果可以通过MAPQ值来提取唯一比对reads，可以用picard、sambamba等软件去除p，最终得到唯一比对且去重复的bam文件。
比对后得到的bam文件可以转化为bigWig（bw）格式，通过可视化软件进行展示。deeptools软件可以实现bw格式转化和可视化展示。首先需要在linux环境中安装deeptools软件，可以用以下命令实现bam向bw格式的转换：
bamCoverage -b A1.bam -o A1.bw
此外，可以使用deeptools软件展示reads在特定区域的分布，如：
computeMatrix reference-point \ # reference-pioint表示计算一个参照点附近的reads分布，与之相对的是scale-regions，计算一个区域附近的reads分布
--referencePoint TSS \#以输入的bed文件的起始位置作为参照点
-S A1.bw \ #可以是一个或多个bw文件
-R gene.bed \ #基因组位置文件
-b 3000 \ #计算边界为参考点上游3000bp
-a 3000 \ #计算边界为参考点下游3000bp，与-b合起来就是绘制参考点上下游3000bp以内的reads分布
-o A1.matrix.mat.gz \ #输出作图数据名称
#图形绘制
plotHeatmap \
-m new_A1.matrix.mat.gz \ #上一步生成的作图数据
-out A1.pdf \ # 输出图片名称
绘图结果展示：
MACS2能够检测DNA片断的富集区域，是ATAC-seq数据call peak的主流软件。峰检出的原理如下：首先将所有的reads都向3'方向延伸插入片段长度，然后将基因组进行滑窗，计算该窗口的dynamic λ，λ的计算公式为：λlocal = λBG（λBG是指背景区域上的reads数目），然后利用泊松分布模型的公式计算该窗口的显着性P值，最后对每一个窗口的显着性P值进行FDR校正。默认校正后的P值（即qvalue）小于或者等于0.05的区域为peak区域。需要现在linux环境中安装macs2软件，然后执行以下命令：
macs2 callpeak \
-t A1.uni.dep.bam \ #bam文件
-n A1 \ # 输出文件前缀名
--shift -100 \ #extsize的一半乘以-1
--extsize 200 \ #一般是核小体大小
--call-summits #检测峰顶信息
注：以上参数参考文献（Jie Wang，et.al.2018.“ATAC-Seq analysis reveals a widespread decrease of chromatin accessibility in age-related macular degeneration.”Nature Communications）
ATAC分析得到的peak是染色质上的开放区域，这些染色质开放区域常常预示着转录因子的结合，因此对peak区域进行motif分析很有意义。常见的motif分析软件有homer和MEME。以homer软件为例，首先在linux环境中安装homer，然后用以下命令进行motif分析：
findMotifsGenome.pl \
A1_peaks.bed \ #用于进行motif分析的bed文件
genome.fa \ #参考基因组fa文件
A1 \ #输出文件前缀
-size given \ #使用给定的bed区域位置进行分析，如果填-size -100,50则是用给定bed中间位置的上游100bp到下游50bp的区域进行分析
homer分析motif的原理及结果参见：http://homer.ucsd.e/homer/motif/index.html
根据motif与已知转录因子的富集情况可以绘制气泡图，从而可以看到样本与已知转录因子的富集显着性。
差异peak代表着比较组合染色质开放性有差异的位点，ChIP-seq和ATAC-seq都可以用DiffBind进行差异分析。DiffBind通过可以通过bam文件和peak的bed文件计算出peak区域标准化的readcount，可以选择edgeR、DESeq2等模型进行差异分析。
在科研分析中我们往往需要将peak区域与基因联系起来，也就是通过对peak进行注释找到peak相关基因。常见的peak注释软件有ChIPseeker、homer、PeakAnnotator等。以ChIPseeker为例，需要在R中安装ChIPseeker包和GenomicFeatures包，然后就可以进行分析了。
library(ChIPseeker)
library(GenomicFeatures)
txdb<- makeTxDbFromGFF(‘gene.gtf’)#生成txdb对象，如果研究物种没有已知的TxDb,可以用GenomicFeatures中的函数生成
peakfile <-readPeakFile(‘A1_peaks.narrowPeak’)#导入需要注释的peak文件
peakAnno <- annotatePeak(peakfile,tssRegion=c(-2000, 2000), TxDb=txdb)
# 用peak文件和txdb进行peak注释，这里可以通过tssRegion定义TSS区域的区间
对于peak注释的结果，也可以进行可视化展示，如：
p <- plotAnnoPie(peakAnno)
通过注释得到的peak相关基因可以使用goseq、topGO等R包进行GO富集分析，用kobas进行kegg富集分析，也可以使用DAVID在线工具来完成富集分析。可以通过挑选感兴趣的GO term或pathway进一步筛选候选基因。
❷ rna-seq不能使用bwa比对吗
对。
_NA测序并不能直接使用DNA测序常用的BWA、Bowtie等比对软件，这是由于真核生物内含子的存在，导致测到的reads并不与基因组序列完全一致（如下图所示），因此需要使用TophatHISATSTAR等专门为RNA测序设计的软件进行比对。
❸ 重测序（RADseq）做群体遗传分析套路
构建的群体，或自然群体，如各地方品种。
提取DNA后，构建文库，简要步骤如下：
① 限制性内切酶TaqI酶切；
② 连接P1接头；
③ DNA随机打断片断化；
④ 目的片段回收与末端修复；
⑤ 连接P2接头；
⑥ RAD片段富集；
⑦ 上机测序。
参考： Rapid and cost-effective polymorphism identification and genotyping using restriction site associated DNA (RAD) markers
根据识别标签序列得到每个个体的测序reads，使用trimmomatic进行过滤（其他质控软件，如fastqc，multiQC等）
设置过滤参数为：SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:50。 过滤标准：两端质量低于5的碱基进行切除，并以5bp为窗口进行滑动过滤，对平均质量低于20的窗口进行切除。
BWA （其他比对软件如bowtie2/soap2/MAQ等）将过滤后的个体clean reads比对到参考基因组序列上。样本比对率反映的是样本测序数据与参考基因组的相似性，覆盖深度和覆盖度能够直接反映测序数据的均一性与参考序列的同源性。
使用GATK（或samtools+bcftools)Haplotype Caller模块进行变异检测，获得群体变异集文件（VCF 格式）。对变异进行过滤：过滤参数为缺失率小于或等于0.2、杂合率小于或等于0.2、最小等位基因频率（MAF） 大于或等于0.05，最终得到高质量的基因型数据。
在获得高质量的标记数据以后，利用vcftools将vcf文件处理得到plink.ped和plink.map文件（整理为plink软件所需格式）。
使用plink 软件随机选择连锁不平衡（LD）小于0.1，且相邻间隔在300kb以上的SNP位点，最后得到一个包含3420 个SNP位点的标记集，一般是生成.bed文件。
1.祖先成分堆叠图
使用ADMIXTURE对此 SNP位点集（bed文件）进行群体结构分析（Structure），利用交叉验证过程确定确定合适的祖先数或亚群（K值）。若不知道理想的K值，可用ADMIXTURE计算，一般当cross-validation error值最低时所对应的K值为最合适的K值。
考虑到样本所归属的分类单元，即看看哪几个物种聚在一起，对合适的K值利用Structure软件（速度慢，其他软件如frappe，ADMIXTURE也可做群体结构图，并且很快）聚类图，一些R包如hapmap也是可以做群体结构图的。
2.PCA
利用GCTA对SNP数据集进行样本的PCA分析（其他软件如EIGENSOFT中的smartpca）。GCTA可以直接读取.bed , .bim , .fam文件，利用–make-grm 生成个体对之间的遗传关系矩阵，并将GRM的下三角元素保存为二进制文件.grm.id , .grm.bin ， .grm.N.bin。使用 –pca 设置要生成主成分的数目，一般来说就可以刻画出群体结构。这一步会生成 .eigenval 和 .eigenvec 两个文件。.eigenval文件为各主成分可解释遗传信息的比例，.eigenvec文件为每个样本在top4主成分上的分解值。
3.系统发育树
构树的方法有非加权分组平均法（UPGMA，已经很少用）、最小进化法（ME）、邻接法（NJ）、最大简约法（MP）、最大似然法（ML）等。
构树软件如FastTree/MEGA/cluster X/phylip，美化可以用FigTree/ggtree/treeview/GraPhIAn。
NJ法是基于最小进化原理经常被使用的一种算法，它不检验所有可能的拓扑结构，能同时给出拓扑结构和分支长度。
GWAS的群体遗传分析也是包含这三个图，RADseq毕竟是简化基因组，得到的SNP有限，做这种群体分析效果肯定没有GWAS好。
❹ RNA-seq中的基因表达量计算和表达差异分析
原文链接： RNA-seq中的基因表达量计算和表达差异分析-生物知识学习 (biotechknowledgestudy.com)
差异分析的步骤：
1）比对；
2） read count计算；
3） read count的归一化；
4）差异表达分析；
背景知识：
1）比对：
普通比对： BWA，SOAP
开大GAP比对：Tophat（Bowtie2）；
2） Read count(多重比对的问题）：
丢弃
平均分配
利用Unique region估计并重新分配
表达量计算的本质
目标基因表达量相对参照系表达量的数值。
参照的本质：
（ 1）假设样本间参照的信号值应该是相同的；
（ 2）将样本间参照的观测值校正到同一水平；
（ 3）从参照的数值，校正并推算出其他观测量的值。
例如：Qpcr:目标基因表达量（循环数）相对看家基因表达量（循环数）；RNA-seq:目标基因的表达量（测序reads数），相对样本RNA总表达量（总测序量的reads数），这是最常用的标准。
归一化的原因及处理原则：
1）基因长度
2）测序量
3）样本特异性（例如，细胞mRNA总量，污染等）前两者使用普通的RPKM算法就可以良好解决，关键是第三个问题，涉及到不同的算法处理。
RNA-Seq归一化算法的意义：
基因表达量归一化：在高通量测序过程中，样品间在数据总量、基因长度、基因数目、高表达基因分布甚至同一个基因的不同转录本分布上存在差别。因此不能直接比较表达量，必须将数据进行归一化处理。
RNA-seq差异表达分析的一般原则
1）不同样品的基因总表达量相似
2）上调差异表达与下调差异表达整体数量相似（上下调差异平衡）
3）在两组样品中不受处理效应影响的基因， 表达量应该是相近的（差异不显着）。
4）看家基因可作为表达量评价依据（ 待定）
不同的算法比较：
以什么数值来衡量表达量：RPKM、FPKM、TPM
以什么作为参照标准：TMM（edgeR软件）、De seq矫正
RPKM：是Reads Per Kilobase per Million mapped reads的缩写，代表每百万reads中来自于某基因每千碱基长度的reads数。
本质：1）以reads数为计算单位；
2）对基因长度（基因间的比较）和总数据量（样本间的比较）做矫正；
1）由于可变剪切，同一基因有效转录区域长度未必相同（这个一般情况下可以不考虑，了解一下：Cufflinks软件考虑了这个问题）优化策略：外显子或转录本水平的表达量分析。
2） 使用reads数计算基因表达量有轻微误差（这里暂不展开，主要了解一下定义）优化策略：FPKM或 TPM
3） mRNA的总量未必相等。
RPKM的优化：FPKm
F = Fragment，即测序片段数量。这些片段都是从完整的cDNA打碎而来的；
本质：以文库中的片段数量为计算单位在Paired-end测序中，一个fragment就是两条PE reads构成的片段。由于是PE比对，理论上比SE比对更可靠。
T = Transcripts
本质：以转录本的条数为计算单位。使用转录本的条数（或者说：转录本的测序深度），代替reads数，在一定条件下定量更准，尤其样本间表达基因总数差异很大的时候（例如，对照样本有1万个基因表达，另外处理组仅有4000个基因表达）。
mRNA总量未必相等
mRNA总量不等——细胞本身不同
例如：活跃组织vs休眠的组织；癌细胞vs正常细胞
mRNA总量不等——污染
例如：核糖体污染外源RNA污染
解决方法——不同算法比较
其中归一化算法介绍：
1）Total Count（TC）：总reads数矫正
2）Upper Quartile（UQ）：上四分之一分位数（总reads）
矫正
3）Median（Med）；中位数（总reads数）矫正
4）Quantile (Q)：基因芯片软件limma中的校正算法；
5）RPKM：总reads数，但引入了基因长度
6）几何平均数：Deseq软件中的算法；
7）TMM：edgeR软件中的算法；
8）RPKM
逻辑1：不同位置数值的稳定性不同
四分位数quartile:将数据按从小到大排列，并分成四等分，这样得到3个分割点，第一个分割点叫做lowerquartile，第二个叫Media，第三个叫Upper quartile
很显然，极大值具有极大不稳定性，而且可能会显着影
响总体之和（假设，我们之中有个马云，我们的总收入
有什么变化？）
所以，Upper quartile和Median的数值，比总表达量之
和更加稳定，更适合作为参照。
逻辑2：表达量居中的基因的表达量值，其数值应该是相似的。
DESeq与edgeR，默认情况下都使用这一的逻辑校正。（DESeq and edgeR Bioconctor packages）
Deseq：异常高表达的基因，会显着影响细胞中的总mRNA的数量。类似的，如果样本中受到不同程度的外源RNA，如病毒、真菌等的污染，也会显着影响样本总mRNA数，导致RPMK值的误差。对于这样的问题，Deseq尝试对数据进行矫正（矫正因子），使表达量处于中间位置的基因表达量应该是基本相同的（即使用表达量处于中间的基因表达量值作为参照，而减少高表达基因的作用）。
Deseq： 校正因子=样本表达中位数/所有样本表达量中位数：回答了一个关键的问题：Deseq不同差异比较组间，计算得到的表达量值不同。因
为样本在变化，“所有样本表达量的中位数”也在变动。RPKM：总表达量为参照
Deseq：中位数为参照
TMM（edgeR）：与Deseq类似，在去除高表达基因和差异最大的基因后，TMM也是要找到一个加权系数，使剩余的基因在被矫正后差异倍数可能小。TMM的加权系数是基于两两样本比较后推算获得的（也就是两组样本的比较，将产生与这次比较相关的加权系数）。然后将所有基因除以这个加权系数，从而保证大部分表达量居中的基因表达量最相似。
不同RNA-seq表达量归一化算法的区别
Deseq类的校正算法：理论上更加稳定；但不同批次的比较会得到不同的表达量值，不利于进行多处理组/批次数据的统一分析（例如，趋势分析、共表达分析）校正会掩盖一些问题（例如：样本污染）
RPKM类的算法： 容易受异常高表达基因、外源污染等的干扰；但也更容易从结果的异常中，发现潜在问题；得到的表达量值是恒定的，多处理组/批次的数据可以合并分析。折中的方法：使用RPKM类的算法，但需要人工检查数据是否
异常。备注： Deseq软件也可以关闭校正的功能。
实际经验总结
总之：从多方面考虑，RPKM类算法，如果合理使用，依然是最优的。具体问题具体分析：在遇到问题的时候，找到问题的来源，从而给出解决方案（没有完美的流程，只有最佳解决方案）
❺ 那令人困惑的比对工具选择啊~
我可能不适合做科研，因为我总是对一些“ 没有必要 ”的事斤斤计较~
刚开始接触二代测序时，是跟着Jimmy大神从RNA-seq开始入门的。那时候使用的比对工具是HISAT2。当时也没怎么细想为什么用这个工具，你说用那我就用吧，于是HISAT2就成了目前我最喜欢的比对工具
但是，随着需要处理不同的数据，我发现HISAT2不能满足我的要求了，我似乎需要考虑其他工具了，但我无从下手，不知道什么工具更适合我。于是一场不太考虑生物学意义，只思考工具特点的战斗便拉开了帷幕……
虽然每次我都说这是废话超多系列，但我是真心想让你们认真看看这部分的内容。毕竟我能力有限，很难用几句简短的话把我想说的事情搞明白，除非我写文言文……
搜了搜各种帖子，发现大家在对比对工具进行比较时，都喜欢将其分为DNA比对工具（DNA-seq）和RNA比对工具（RNA-seq）。仔细思考你会发现，它们的区别仅在于是否会考虑跨外显子的比对（即：是否会将没有比对上的reads劈开，对劈开后的两部分再次比对）。
随着现在各种seq的出现，我们已经不能简单的根据是比对DNA还是RNA来判断工具的选择，而是要判断reads的比对是否需跨外显子。比如PRO-seq/GRO-seq，它们在建库时捕获的RNA，但是它们并不需要考虑跨外显子的比对。
鬼扯了这么多，简单总结一下各种类型的常用工具都有哪些：
bowtie出现在上古时期（就是很久远的意思了），那个时候测序行业的发展还不成熟，序列长度普遍在50bp以下，因此bowtie的出现就是为了满足长度在50bp以下的reads的比对。官方称其可以把短的DNA序列（35bp）快速的比对到人类基因组上。
而bowtie2的出现则弥补了bowtie的短板，bowtie2擅长比对50-100bp长的reads，长度甚至长达1000。它适合比对那些比较长的基因组，如哺乳动物基因组。
结论：bowtie和bowtie2，是两个不同类型的比对工具，bowtie2并非是bowtie的升级。尺有所长寸有所短，bowtie适合长度在50b长度以内的reads比对，而bowtie2适合50-100b，甚至更长的reads比对。但是这两个都属DNA-seq比对工具
Tophat/Tophat2工具本身不能进行比对，它是通过调用bowtie/bowtie2进行比对的。划重点， bowtie2不是bowtie的升级版，但是Tophat2是Tophat2的升级版 。因此Tophat只可以调用bowtie，而Tophat2不仅可以调用bowtie2（默认）还可以更改设置调用bowtie。
Tophat/Tophat2调用bowtie/bowtie2后，会首先使用bowtie/bowtie2对序列进行比对，对于那些没有比对上的，会考虑其跨外显子的可能性，将reads劈开重新比对。
如果你去bowtie/bowtie2/Tophat的官网仔细观察，你会发现，bowtie和bowtie2各自有自己的官网，有专属于自己的介绍。而Tohat就不同了，它只有一个，仅仅是在2012年4月9日Tophat发布了2.0.0版本，宣布支持bowtie2的比对。而我们通常也将支持bowtie2版本的Tophat称之为Tophat2。
此外，如果你够无聊，你在它们的主页上下扒拉扒拉，你会发现无论是bowtie还是bowtie2在2019年仍然是有更新的。但是Tophat到了2016年2月便停止了更新……这是为什么呢？请继续往下看
Tophat2的原作者们也不知道是出于什么考虑，不再更新Tophat2，转而开发了一个新的比对工具HISAT2，更是推荐人们使用HISAT2，声称其速度更快，内存占用率更小，准确率更高。
此外，HISAT2不仅支持RNA-seq的比对还支持DNA-seq比对，唯一需要做的就是加上一个参数 --no-spliced-alignment 。但是就目前来看，大部分人都是使用HISAT2做RNA-seq，没人使用它做DNA-seq
其实没有太多需要说的，是我最早知道的比对工具，属于DNA-seq比对工具。大部分搞全基因组或者全外的似乎都使用BWA作比对。当时本科时期有门课叫做计算生物学，当时学习BWT算法，还经常把算法名称和比对工具名称搞混。
值得一提的是这个工具是李恒开发的，如果你不知道李恒，你最起码也应该知道SAMtools工具，这东西也是李恒开发的。一个强的可怕的男人……
最早听到这个工具的名称是在研一时期，的确有点孤陋寡闻。当时以为是一个很小众的野鸡工具，后来发现身边蛮多人用这个工具的，于是就带着好奇心上网搜索它的资料。
不搜不知道，一搜吓一跳。啧啧啧，ENCODE皇家御用的RNA-seq比对工具，真香……
可能是因为RNA-seq分析比较大众，因此大部分的比对工具都是利用RNA-seq的效率进行比较。
《Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis》这篇文章可以称之为史上最全RNA-seq测评。对于RNA-seq的方方面面都进行了比较。因为太全面，反而有点杂乱，因此我只关注了我感兴趣的一些地方。
无论是HISAT2还是STAR，对于Tophat来说都有很大的优势，更何况Toph还不再继续更新，这就更给力我们不再使用它的理由。至于HISAT2，他的junction正确率最高，但是灵敏度相对较低。而STAR灵敏度更高，但是会有许多包含soft-clip的低质量比对。此外，STAR的unique mapping比例最高，它对于双端测序的reads，要么全部比对上，要么全部抛弃，不会像像TopHat和HISAT2一样只比对上某一个reads
最后这篇文章还给了一个它认为比较好的RN-seq组合方式：
前一阵的龙星课程针对RNA-seq给出了另一个组合方式：SATR+RSEM，其中STAR既可以比对也可以用于定量（count）
具体选用哪一个组合，看习惯，看眼缘，看心情……
❻ 如何使用bwa只输出一条比对结果
首先要把所有的序列复制到windows自带的记事本中,全部以fasta格式存到同一个文件中,保存成*.txt,序列内部最好不要有空格或者换行.
序列格式举例如下：
“>序列1
ATCG.ATCG
>序列2
ATCG.ATCG
>序列3
ATCG.ATCG
>序列4
ATCG.ATCG”
然后用Bioedit打开刚才保存的文件*.txt,点击窗口上方的Accessory application菜单,再点击clustelW multiple alignment,这时候会弹出一个窗口,直接选择Run clustalW,又弹出一个窗口,选择ok,等待结果就行了.
最后的比对结果可以保存成*.fas格式,适用于各种分析.