Biostar第七课02 align到基因组
这个是重点中的重点,虽然是最基础的,但是这个的结果将直接影响下游的多种分析,可以说是核心文件。
怎样去选择参数,这个貌似要看自己的观点,其实就是告诉软件怎样打分,或者更确切的说,mismatch扣几分,deltion扣几分,开个gap扣几分
然后软件根据你制定的规则来找出得分最高的align的位置
几乎所有的alinger的默认选择都是EDNAFULL标准
这个标准可以这样下载
curl -O ftp://ftp.ncbi.nlm.nih.gov/blast/matrices/NUC.4.4
#打分矩阵
A T G C
A 5 -4 -4 -4
T -4 5 -4 -4
G -4 -4 5 -4
C -4 -4 -4 5
关于怎样选择合适的矩阵,有这个WIlliam Pearson 牛逼哄哄的论文可以参考
这里还有各种矩阵可供选择
不过得小心一点,有一些标准化过了,有一些没有标准化,虽然名字可能看起来都差不多
Gap penalty
这个也是比较重要的参数,penalty的标准有很多,一般来说软件都用Affine的gap penalty,这个考虑两个参数。第一个就是gap的open,也就是发现一个gap。第二个是gap的extension,也就是这个gap的长度。如果你觉得,出现gap比较讨厌,就可以将open的罚分设的高一点。如果觉得出现比较长的gap很讨厌,就把extension的罚分搞高一点就好。但是一般来说,open的罚分会比较高,extension的罚分相对来说小很多。
全局对比
目的是看两条序列的相似程度
这里两条序列中每一个碱基都得对齐到另一条序列的相同碱基,或者是错配碱基上,或者是gap上。划重点,每一条每一个碱基
工具是EMboss套装,工具的开发者本意是帮助生信不熟的生物僧来使用,所以用起来会有很多的提示信息,可能会比较烦人,所以加个参数‘-filter’, 世界顿时就清净了许多。
里面的全局对比的命令使用gapopen罚分是10分一个,软件觉得开一个gap就丢很多分,所以就倾向于少开gap,宁可准确对齐的碱基少一点,也比开个gap罚分小
但是当我们用 -gapopen 7 将罚分降低一点,软件顿时就舒服了很多,开gap也不那么谨慎了,所以对齐的碱基可能更多一些,但是相应的gap也多了很多。
还可以用‘-data’ 指定打分的矩阵文件
一般来说软件可能会采用 free-end-gap,就是说最末尾有gap罚分不那么高,
局部对比
主要是用来寻找两条序列中相似度最大的小区域,划重点,其中的一小段区域
算法主要是寻找两条序列中各自的一小段序列,这两小段序列相似度特别高
打分矩阵不同,结果很不一样
# DNA矩阵
# This matrix is the "standard" EDNAFULL substitution matrix.
wget ftp://ftp.ncbi.nlm.nih.gov/blast/matrices/NUC.4.4
#蛋白矩阵
# Get the BLOSUM30, BLOSUM62, and BLOSUM90 matrices.
wget ftp://ftp.ncbi.nlm.nih.gov/blast/matrices/BLOSUM30
wget ftp://ftp.ncbi.nlm.nih.gov/blast/matrices/BLOSUM62
wget ftp://ftp.ncbi.nlm.nih.gov/blast/matrices/BLOSUM90
半全局对比
这个是我们测序分析的时候最常用的算法,主要用来将一小段短序列对比到很长的序列中,也就是将测序的reads对比到参考基因组中
准确配对,配对误差,软件的局限性
看了这么多,其实最关心的就是,到底准不准,到底结果能不能还原出真实的生物学事件。
作者做了个小实验,在原始序列里面的长串的C两边各插入一个C,结果全局对比到原始序列时,第一位置错误,第二将两个C的insertion报告为一个CT插入后紧跟一个T错配。
原因就是这一串串的C里面给出的信息太少,aligner从中获得不到任何的纠错区分能力,aligner只能寻求数学上的最优解,而不能找出更加符合生物学解释的方案。
这是不是说数学上的最优解不一定是生物学的最优解,所以生物僧做好实验还是意义重大的。
以上的问题可以用一个参数很轻易的解决,就是gapopen=9,也就是一个gap不要罚十分,罚九分就好。
但这个绝对不是说以后就要用gapopen=9 而不是10这个默认参数,调整参数的影响是很深远的,直接会影响到本来会正确的align的reads。使用9这个参数的结果就是,两个正确的碱基配对得分5+5=10分,而一个gap只罚分9分,小于10,所以align在遍历基因组的过程中,只要找到一个gap跟着两个碱基的match就会得分正值,判定为正确配对。在某些特定的情况下,可能挺好,但是绝大多数情况不太行。
以上的错误,并不是由于参数设置错误所导致的,问题的本质在于连续多个C导致的均聚区域给出配对信息太少,软件无法根据这么少的信息纠正错误。
配对可靠性取决于序列本身自己的复杂程度以及能够给出的信息量
复杂度比较低的区域通常给出的align的结果不是很可靠,需要额外的验证
多条序列align
mafft --clustalout small.fa > alignment.maf
KR105345.1 -------ataattttcctctcattgaaatttatatcggaatttaaattgaaattgttact
KR105328.1 --gattaataattttcctctcattgaaatttatatcggaatttaaattgaaattgttact
KR105323.1 --gattaataattttcctctcattgaaatttatatcggaatttaaattgaaattgttact
KR105302.1 --gattaataattttcctctcattgaaatttatatcggaatttaaattgaaattgttact
KR105295.1 ---attaataattttcctctcattgaaatttatatcggaatttaaattgaaattgttact
KR105294.1 --gattaataattttcctctcattgaaatttatatcggaatttaaattgaaattgttact
KR105282.1 --gattaataattttcctctcattgaaatttatatcggaatttaaattgaaattgttact
KR105266.1 ---attaataattttcctctcattgaaatttatatcggaatttaaattgaaattgttact
KR105263.1 aagattaataattttcctctcattgaaatttatatcggaatttaaattgaaattgttact
KR105253.1 ---attattaatyttcctctcattgaaatttatatcggaatttaaattgaaattgttact
**** ***********************************************
还可以用clustal-omega 来多序列对比