family_pca.R

set.seed(42)

# Num individuals
par( mfrow=c(1,4) )
# Num markers
M = 500
# Minor allele frequency
maf = 0.5

LABEL = TRUE

title_prefix = c("A","B","C")
all_N = c(1,5,50)

for ( i in 1:3 ) {
  N = all_N[i]
  mate1_hap_a = matrix( rbinom(N*M,1,maf) , nrow=N , ncol=M )
  mate1_hap_b = matrix( rbinom(N*M,1,maf) , nrow=N , ncol=M )
  
  mate2_hap_a = matrix( rbinom(N*M,1,maf) , nrow=N , ncol=M )
  mate2_hap_b = matrix( rbinom(N*M,1,maf) , nrow=N , ncol=M )
  
  n_children = 5
  all_children_hap_a = matrix( NA , nrow=n_children , ncol=M )
  all_children_hap_b = matrix( NA , nrow=n_children , ncol=M )
  
  for ( n in 1:n_children ) {
    # generate offspring
    child_hap_a = mate1_hap_a
    sampled_vars_mate1 = matrix( as.logical(rbinom(N*M,1,0.5)) , nrow=N , ncol=M )
    child_hap_a[ !sampled_vars_mate1 ] = mate1_hap_b[ !sampled_vars_mate1 ]
    
    child_hap_b = mate2_hap_a
    sampled_vars_mate2 = matrix( as.logical(rbinom(N*M,1,0.5)) , nrow=N , ncol=M )
    child_hap_b[ !sampled_vars_mate2 ] = mate2_hap_b[ !sampled_vars_mate2 ]
    
    all_children_hap_a[n,] = child_hap_a[1,]
    all_children_hap_b[n,] = child_hap_b[1,]
  }
  
  all_genos = all_children_hap_a + all_children_hap_b
  all_genos = rbind( all_genos , mate1_hap_a + mate1_hap_b , mate2_hap_a + mate2_hap_b )
  
  # drop monomorphic
  keep = apply(all_genos,2,sd) != 0
  all_genos = all_genos[ , keep ]
  
  pc <- prcomp( t( scale(all_genos) ) , center = FALSE, scale. = FALSE )
  varexp = pc$sdev^2 / sum(pc$sdev^2)
  
  if ( LABEL ) {
    point_clr = c( rep("black",n_children), "black" , rep("#31a354",(N-1)) , "black" , rep("#31a354",(N-1))  )
    point_typ = c( rep(19,n_children),1,rep(19,N-1),1,rep(19,N-1))
    plot( pc$rotation[,"PC1"] , pc$rotation[,"PC2"] , col = point_clr  , las=1 , cex=0.75 , pch=point_typ , xlab="PC1" , ylab="PC2" , main=paste(title_prefix[i],": ",N*2," unrelated samples",sep='') )
    
  } else {
    point_clr = "black"
    point_typ = 19
    plot( pc$rotation[,"PC1"] , pc$rotation[,"PC2"] , col = point_clr  , las=1 , cex=0.75 , pch=point_typ , xlab="PC1" , ylab="PC2" , main=title_prefix[i] )
    
  }

  if ( N == 1 && LABEL ) {
    legend("topleft",legend=c("unrel parents","children","unrel others") , bty="n" , cex=1 , pch=c(1,19,19) , col=c("black","black","#31a354"))
  }
}

pc <- prcomp( t( scale(all_genos[ ((n_children+1):nrow(all_genos)) ,]) ) , center = FALSE, scale. = FALSE )
varexp = pc$sdev^2 / sum(pc$sdev^2)

if ( LABEL ) {
  point_clr = c( "black" , rep("#31a354",(N-1)) , "black" , rep("#31a354",(N-1))  )
  point_typ = c( 1,rep(19,N-1),1,rep(19,N-1))
  plot( pc$rotation[,"PC1"] , pc$rotation[,"PC2"] , col = point_clr  , las=1 , cex=0.75 , pch=point_typ , xlab="PC1" , ylab="PC2" , main="D: Only unrelated" )
  
} else {
  point_clr = "black"
  point_typ = 19
  plot( pc$rotation[,"PC1"] , pc$rotation[,"PC2"] , col = point_clr  , las=1 , cex=0.75 , pch=point_typ , xlab="PC1" , ylab="PC2" , main="D" )
  
}