BSA分析中MutMap 分析:95% 99%阈值置信区间是如何计算的

BSA分析中MutMap 分析:95% 99%阈值置信区间是如何计算的


R语言代码有兴趣可以了解一下:

###########genotype#############################
genotype<-function(){
count<-0


if (population_structure=="RIL"){
x<-runif(1) 
if (x<=0.5){
count<-1
}else{
count<-0
}

}else{
for(i in 1:2){
x<-runif(1) 
if (x<=0.5){
number<-0.5
}else{
number<-0
}
if(number == 0.5){
count<- count+0.5
}
}
}
return(count)
}
############################################################

###########caluclate of genotype ratio#########################
individuals_genotype<-function(number_of_total_individuals){

ratio_of_genotype<-c()
for(i in 1:number_of_total_individuals){
ratio_of_genotype<-c(ratio_of_genotype,genotype())

}
return(mean(ratio_of_genotype))
}
############################################################

###########SNP_index_caluclation#########################
snp_index<-function(read_depth,ratio_of_genotype_in_the_population){
x1<-rbinom(1,read_depth,ratio_of_genotype_in_the_population)
return(x1/read_depth)
}
############################################################
####################################

####################################
Arg<-commandArgs(TRUE)
###########input###############################################
population_structure<-"F2"
individual_analysis<- c(Arg[1])
reprication<-c(Arg[2])
filter_value<-c(Arg[3])
depth_analysis<-c(1:300)
###########input###############################################



for (key_individual in individual_analysis){
   individual_number<-key_individual
    
    depth_data<-c()
    p_h_data_95<-c()
    p_h_data_99<-c()
    for (key_depth in depth_analysis){
            depth_data<-c(depth_data,key_depth)
            depth<-key_depth
        sum_snp_index<-c()    
        for(i in 1:reprication){    
        ##########gene_frequency######################
ratio_of_genotype_in_the_population<-individuals_genotype(key_individual)
a_snp_index<-snp_index(key_depth,ratio_of_genotype_in_the_population)
if(a_snp_index >= filter_value){
sum_snp_index<-c(sum_snp_index,a_snp_index)
}
##########gene_frequency######################
        }
        order_sum_snp_index<-sort(sum_snp_index)
        length_sum_snp_index<-length(sum_snp_index)
        ##########snp_index_probabirity_0.05######################       
snp_cutoff_up_0.95<-order_sum_snp_index[ceiling(0.95*length_sum_snp_index)]
        p_h_data_95<-c(p_h_data_95,snp_cutoff_up_0.95)   
        ##########snp_index_probabirity_0.05######################
        
        ##########snp_index_probabirity_0.01###################### 
        
        if (ceiling(0.99*length_sum_snp_index)<length_sum_snp_index){
snp_cutoff_up_0.99<-order_sum_snp_index[ceiling(0.99*length_sum_snp_index)]
}else{
snp_cutoff_up_0.99<-order_sum_snp_index[length_sum_snp_index]
}
        p_h_data_99<-c(p_h_data_99,snp_cutoff_up_0.99) 
        ##########snp_index_probabirity_0.01######################      
               
            
     
     
    }
    
        
    FINAL_DATA<-data.frame(DEPTH=depth_data,P_H_95=p_h_data_95,P_H_99=p_h_data_99)
    
    table_name<-paste("./",individual_number,sep="")
    table_name<-paste(table_name,"individuals.txt",sep="")
    write.table(FINAL_DATA,table_name,sep="\t", quote=F, append=F,row.name=F)
    
}  
    
    
    
    



  • 发表于 2020-07-21 13:02
  • 阅读 ( 108 )
  • 分类:重测序

0 条评论

请先 登录 后评论
omicsgene
omicsgene

生物信息

327 篇文章

作家榜 »

  1. omicsgene 327 文章
  2. 安生水 217 文章
  3. Daitoue 167 文章
  4. 生物女学霸 120 文章
  5. CORNERSTONE 72 文章
  6. 生信老顽童 48 文章
  7. 红橙子 47 文章
  8. landy 37 文章