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

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

这里R语言源代码,有兴趣的可以研究看看



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


###########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_in_A){
x1<-rbinom(1,read_depth,ratio_of_genotype_in_the_population_in_A)
return(x1/read_depth)
}
############################################################
####################################





for (key_individual in individual_analysis){
   individual_number<-key_individual
    
    depth_data<-c()
    p_l_data_95<-c()
    p_h_data_95<-c()
    p_l_data_99<-c()
    p_h_data_99<-c()
    for (key_depth in depth_analysis){
        depth_data<-c(depth_data,key_depth)
        depth<-key_depth
    
        data_of_delta_snp_index<-c()
        
        for(i in 1:reprication){    
        ##########gene_frequency######################
ratio_of_genotype_in_the_population_in_A<-individuals_genotype(key_individual)
Snp_index_of_A<-snp_index(key_depth,ratio_of_genotype_in_the_population_in_A)

ratio_of_genotype_in_the_population_in_B<-individuals_genotype(key_individual)
Snp_index_of_B<-snp_index(key_depth,ratio_of_genotype_in_the_population_in_B)

if(Snp_index_of_A >= filter_value | Snp_index_of_B >=filter_value){
delta_snp_index<-Snp_index_of_A-Snp_index_of_B
data_of_delta_snp_index<-c(data_of_delta_snp_index,delta_snp_index)
}

##########gene_frequency######################
        }
        
        order_data_of_delta_snp_index<-sort(data_of_delta_snp_index)
        length_data_of_delta_snp_index<-length(data_of_delta_snp_index)
        ##########snp_index_probabirity_0.05######################       
        snp_cutoff_low_0.025<-order_data_of_delta_snp_index[floor(0.025*length_data_of_delta_snp_index)]
snp_cutoff_up_0.975<-order_data_of_delta_snp_index[ceiling(0.975*length_data_of_delta_snp_index)]
p_l_data_95<-c(p_l_data_95,snp_cutoff_low_0.025)
        p_h_data_95<-c(p_h_data_95,snp_cutoff_up_0.975)   
        ##########snp_index_probabirity_0.05######################
        
        ##########snp_index_probabirity_0.01###################### 
        if (floor(0.005*length_data_of_delta_snp_index)>0){
        snp_cutoff_low_0.005<-order_data_of_delta_snp_index[floor(0.005*length_data_of_delta_snp_index)]
        }else{
        snp_cutoff_low_0.005<-order_data_of_delta_snp_index[1]
        }      
        
        if (ceiling(0.995*length_data_of_delta_snp_index)<length_data_of_delta_snp_index){
snp_cutoff_up_0.995<-order_data_of_delta_snp_index[ceiling(0.995*length_data_of_delta_snp_index)]
}else{
snp_cutoff_up_0.995<-order_data_of_delta_snp_index[length_data_of_delta_snp_index]
}
p_l_data_99<-c(p_l_data_99,snp_cutoff_low_0.005)
        p_h_data_99<-c(p_h_data_99,snp_cutoff_up_0.995) 
        ##########snp_index_probabirity_0.01######################      
               
            
     
     
    }
            
    FINAL_DATA<-data.frame(DEPTH=depth_data,P_L_95=p_l_data_95,P_H_95=p_h_data_95,P_L_99=p_l_data_99,P_H_99=p_h_data_99)
    table_name<-paste("./",population_structure,"_",individual_number,"_individuals.txt",sep="")
    
    write.table(FINAL_DATA,table_name,sep="\t", quote=F, append=F,row.name=F)
    
}  
    

    

    

    模拟好的数据:simulation_result.zip





  • 发表于 2020-07-21 12:51
  • 阅读 ( 4824 )
  • 分类:重测序

2 条评论

请先 登录 后评论
omicsgene
omicsgene

生物信息

698 篇文章

作家榜 »

  1. omicsgene 698 文章
  2. 安生水 343 文章
  3. Daitoue 167 文章
  4. 生物女学霸 120 文章
  5. xun 78 文章
  6. 红橙子 78 文章
  7. rzx 74 文章
  8. CORNERSTONE 72 文章