有时候我们要利用C++处理bam文件,比如取出read1要么read2等适合特定条件的班,根据cigar值对班指定位置的碱基进行统计要对班进行拍卖并出口相当于,这时我们可使用htslib库。htslib可以据此来处理SAM,
BAM,CRAM 和VCF文件,是samtools、bcftools的核心库。

实例

  生物数据解析中,我们想查PCR扩增出来的扩增子的测序深度以内的出入,但不同的扩增子的扩增效率受到GC含量之影响,因此我们首先应当解除掉GC含量对扩增子深度的震慑。

参考资料

htslib
sam.h文件:https://github.com/samtools/htslib/blob/develop/htslib/sam.h
htslib
sam文件格式说明:https://github.com/samtools/hts-specs/blob/master/SAMv1.pdf

R语言代码

 loess(formula, data, weights, subset, na.action, model = FALSE,
       span = 0.75, enp.target, degree = 2,
       parametric = FALSE, drop.square = FALSE, normalize = TRUE,
       family = c("gaussian", "symmetric"),
       method = c("loess", "model.frame"),
       control = loess.control(...), ...)

  formula是公式,比如y~x,可以输入1交4只变量;
  data是推广正变量的数据框,如果data为空,则于环境面临找找;
  na.action指定对NA数据的拍卖,默认是getOption(“na.action”);
  model是否返模型框;
  span是alpha参数,可以决定平滑度,相当给面所陈述之f,对于alpha小于1的上,区间涵盖alpha的触及,加权函数为立方加权,大于1时,使用有的接触,最深距为alpha^(1/p),p
为讲变量;
  anp.target,定义span的准备方式;
  normalize,对多变量normalize到同一scale;
  family,如果是gaussian则使最小二乘机法,如果是symmetric则使用对姑函数进行再下降的M估计;
  method,是适应模型或仅提取模型框架;
  control进一步再尖端的主宰,使用loess.control的参数;
  其它参数请自己参见manual并且查找资料

loess.control(surface = c("interpolate", "direct"),
          statistics = c("approximate", "exact"),
          trace.hat = c("exact", "approximate"),
          cell = 0.2, iterations = 4, ...)

  surface,拟合表面是于kd数进行插值还是进行规范计算;
  statistics,统计数据是纯正计算还是近似,精确计算好缓慢
  trace.hat,要跟的平整的矩阵精确计算还是近似?建议采取过1000个数据点逼近,
  cell,如果通过kd树最可怜的触及开展插值的好像。大于cell
floor(nspancell)的点让划分。
  robust fitting使用的迭代次数。

predict(object, newdata = NULL, se = FALSE,
    na.action = na.pass, ...)

  object,使用loess拟合出来的目标;
  newdata,可选取数据框,在内部找变量并拓展展望;
  se,是否算标准误差;
  对NA值的处理

cigar值存储形式

32位int,通过bam_get_cigar获得地方,aln->core.n_cigar返回cigar
operation的个数

  • 低 4位存储cigar operation;通过函数bam_cigar_op()获得operation
    统计 1
  • 赛28各类存储cigar值的长短;通过函数,bam_cigar_oplen获得
LOESS平滑方法

  1.
以x0也中心确定一个距离,区间的宽度可以活掌握。具体来说,区间的宽窄在q=fn。其中q是介入有回归观察值的个数,f是在座一些回归观察值的个数占观察值个数的比例,n是观察值的个数。在实际上利用中,往往先选定f值,再根据f和n确定q的取值,一般景象下f的取值在1/3交2/3内。q与f的取值一般从不规定的守则。增大q值或f值,会导致平滑值平滑程度大增,对于数据中前在的轻变化模式则分辨率低,但噪声小,而针对性数码中大的浮动模式之展现虽然较好;小之q值或f值,曲线粗糙,分辨率高,但噪声大。没有一个专业的f值,比较明智之做法是连的调节比较。
  2.
概念区间内所有点的权数,权数由权数函数来规定,比如立方加权函数weight =
(1 –
(dist/maxdist)^3)^3),dist为距离x的偏离,maxdist为距离内距离x的绝可怜距。任一点(x0,y0)的权数是权数函数曲线的冲天。权数函数应包括以下三独面特征:(1)加权函数上之点(x0,y0)具有最老权数。(2)当x离开x0(时,权数逐渐滑坡。(3)加权函数以x0为中心对称。
  3.
对准区间内的散点拟合一修曲线y=f(x)。拟合的直线反映直线关系,接近x0的触及于直线的拟合中自至重要的来意,区间外之接触它的权数为零星。
  4.
x0的平滑点就是x0在拟合出来的直线上的拟合点(y0,f(
x0))。
  5. 针对性持有的点请求来平滑点,将平滑点连接就收获Loess回归曲线。

seq存储形式

8各int,4各存储一个碱基,1,2,4,8,15分级代表A、C、G、T、N,高四位存储坐标数较小之碱基,可由此bam_seqi(seq,i)遍历。

  当我们怀念研究不同sample的有变量A之间的区别经常,往往会盖任何一些变量B对该变量的原来影响,而影响不同sample变量A的比较,这个时需要对sample变量A进行规范后才能够开展较。标准化的方法是针对sample
的 A变量和B变量进行loess回归,拟合变量A关于变量B的函数
f(b),f(b)则意味着于B的熏陶下A的驳斥取值,A-f(B)(A对f(b)残差)就得错过掉B变量对A变量的震慑,此时残差值就可当条件的A值在不同sample之间展开较。

#include <stdio.h>
#include <stdlib.h>
#include <htslib/sam.h>

using namespace std; 

#define bam_is_read1(b) (((b)->core.flag&BAM_FREAD1) != 0)

uint8_t Base[16] = {0,65,67,0,71,0,0,0,84,0,0,0,0,0,0,78};

int main(int argc, char **argv)
{
    bam_hdr_t *header;
    bam1_t *aln = bam_init1();

    samFile *in = sam_open(argv[1], "r");
    htsFile *outR1 = hts_open(argv[2], "wb");
    header = sam_hdr_read(in);
    if (sam_hdr_write(outR1, header) < 0) {
        fprintf(stderr, "Error writing output.\n");
        exit(-1);
    }
    uint8_t *seq;
    int32_t lseq;
    uint32_t *cigar;
    char* qname;
    while (sam_read1(in, header, aln) >= 0) {
        if (bam_is_read1(aln)){
            sam_write1(outR1, header, aln);
        }
        else {
            seq = bam_get_seq(aln);
            lseq = aln->core.l_qseq;
            qname  = bam_get_qname(aln);
            printf("%s\n",qname);
            cigar = bam_get_cigar(aln);
            for(int i=0; i < aln->core.n_cigar;++i){
                int icigar = cigar[i];
                printf("%d%d\n",bam_cigar_op(icigar),bam_cigar_oplen(icigar));
            }
            for(int i=0; i < lseq;++i){
               printf("%c", Base[bam_seqi(seq, i)]);
            }
            printf("\n");

        }
    }
    sam_close(in);
    sam_close(outR1);
}

数据

amplicon
测序数据,处理后拿走的每个amplicon的纵深,每个amplicon的GC含量,每个amplicon的尺寸
统计 2
优先用loess进行曲线的拟合

gcCount.loess <- loess(log(RC+0.01)~GC,data=RC_DT,control = loess.control(surface = "direct"),degree=2)

绘画出拟合出来的曲线

predictions1<- predict (gcCount.loess,RC_DT$GC)
#plot scatter and line 
plot(RC_DT$GC,log(RC_DT$RC+0.01),cex=0.1,xlab="GC Content",ylab=expression(paste("log(NRC"["lib"],"+0.01)",sep="")))
lines(RC_DT$GC,predictions1,col = "red")

统计 3

取残差,去除GC含量对纵深的熏陶

#sustract the influence of GC
resi <- log(RC_DT$RC+0.01)-predictions1
RC_DT$RC <- resi
setkey(RC_DT,GC)

此时RC_DT$RC就是normalize之后的RC
画画图展示nomalize之后的RC,并拿拟合的loess曲线和normalize之后的数额保存

#plot scatter and line using Norm GC data
plot(RC_DT$GC,RC_DT$RC,cex=0.1,xlab="GC Content",ylab=expression("NRC"["GC"]))
gcCount.loess <- loess(RC~GC,data=RC_DT,control = loess.control(surface = "direct"),degree=2)
save(gcCount.loess,file="/home/ywliao/project/Gengyan/gcCount.loess.Robject")
predictions2 <- predict(gcCount.loess,RC_DT$GC)
lines(RC_DT$GC,predictions2,col="red")
save(RC_DT,file="/home/ywliao/project/Gengyan/RC_DT.Rdata")

统计 4

理所当然,也想看一下amplicon 长度len 对RC的震慑,不过影响不很
统计 5

全方位代码如下(经过改,可能跟方了配合):

library(data.table)

load("/home/ywliao/project/Gengyan/RC_DT.Rdata")
RRC_DT <- RC_DT[Type=="WBC" & !is.na(RC),]

lst <- list()
for (Samp in unique(RC_DT$Sample)){
RC_DT <- RRC_DT[Sample==Samp]
####loess GC vs RC####
gcCount.loess <- loess(log(RC+0.01)~GC,data=RC_DT,control = loess.control(surface = "direct"),degree=2)
predictions1<- predict (gcCount.loess,RC_DT$GC)
#plot scatter and line 
#plot(RC_DT$GC,log(RC_DT$RC+0.01),cex=0.1,xlab="GC Content",ylab=expression(paste("log(NRC"["lib"],"+0.01)",sep="")))
#lines(RC_DT$GC,predictions1,col = "red")
#sustract the influence of GC
resi <- log(RC_DT$RC+0.01)-predictions1
RC_DT$NRC <- resi
setkey(RC_DT,GC)
#plot scatter and line using Norm GC data
#plot(RC_DT$GC,RC_DT$NRC,cex=0.1,xlab="GC Content",ylab=expression("NRC"["GC"]))
gcCount.loess <- loess(NRC~GC,data=RC_DT,control = loess.control(surface = "direct"),degree=2)
predictions2 <- predict(gcCount.loess,RC_DT$GC)
#lines(RC_DT$GC,predictions2,col="red")
lst[[Samp]] <- RC_DT
}
NRC_DT <- rbindlist(lst)
save(RC_DT,file="/home/ywliao/project/Gengyan/NRC_DT.Rdata")

####loess len vs RC###
setkey(RC_DT,Len)
len.loess <- loess(RC_DT$NRC~RC_DT$Len, control = loess.control(surface = "direct"),degree=2)
predictions2<- predict (len.loess,RC_DT$Len)
#plot scatter and line 
plot(RC_DT$Len,RC_DT$NRC,cex=0.1,xlab="Length",ylab=expression(paste("log(RC"["GC"],"+0.01)",sep="")))
lines(RC_DT$Len,predictions2,col = "red")

Loess局部加权多项式回归

  LOWESS最初由Cleveland
提出,后而受Cleveland&Devlin及其他多人数进步。在R中loess
函数是坐lowess函数为根基的又复杂功能再强的函数。主要考虑为:在数据集合的各个一点用低维多项式拟合数据点的一个子集,并估计该点附近自变量数据点所对应之以变量值,该多项式是为此加权最小二趁法来拟合;离该点越远,权重越小,该点的回归函数值就是这局部多项式来抱,而用于加权最小二随着回归的数据子集是出于多年来邻近方法确定。
  最为要命亮点:不待先设定一个函数来对拥有数据拟合一个模。并且可以对同一数据进行反复例外的拟合,先对某个变量进行拟合,再指向任何一样变量进行拟合,以探索数据被或者在的某种关联,这是平常的回归拟合无法到位的。