<commit_before># Analysis of Genetic Data 1:<br>Making predictions using PCA

Here we further explore principal components by computing the
"projection" of a new set of genotype samples (a "test set") onto the
PCs computed previously using a different set samples (a "training set").

To be able to work with new genotype samples, we need to make sure
that they are represented, or encoded, in the genotype matrix *in the
exact same way.* In particular, we need to make sure that a "0"
genotype in the test samples means the same thing as a "0" genotype in
the training samples, and the same for "1" and "2" genotypes. We can
enforce the same encoding this using the `--recode-allele` flag in
PLINK.

```bash
module load plink/1.90
cut -f 2,5 1kg_train.bim > alleles.txt
plink --file 1kg_test --recode A-transpose spacex \
  --recode-allele alleles.txt --out 1kg_test
```

We now have the genotype data stored as a matrix with entries that are
0, 1, 2 or missing (NA). Next, we open up R again,

```bash
cd ~/git/genetic-data-analysis-rcc-1/code
R --no-save
```

and in R we enter the following:

```R
# Load the functions for reading data from various files.
source("read.data.R")

# Read the genotypes of the test samples from the .traw file. 
X <- read.geno.traw(fam.file  = "../data/1kg_test.fam",
                    geno.file = "../data/1kg_test.traw")
print(dim(X))

# Read the PCA rotation generated by MATLAB script geno_pca.m.
R <- read.rot.matrix("../results/1kg_train_rot.txt")
print(dim(R))

# Read the mean genotypes outputed by geno_pca.m.
y <- read.mean.genotypes("../results/1kg_train_mean.txt")
```

The genotype matrix we just loaded should have one row for each of the
29 test samples and one column for each of the 156,923 SNPs (same as
the training samples). The rotation matrix is a p x k matrix, where k
is the number of PCs computed.

A simple matrix multiplication completes the projection of the test
samples onto the embedding defined by the PCs.

```R
source("misc.R")
X           <- center.cols(X,y)
X[is.na(X)] <- 0
pc.test     <- X %*% R
```

The only annoying problem was the small number (<1%) of missing
genotypes that make the matrix multiplication less straightforward.
Here we simply set the missing entries to zero after centering the
columns of the matrix. Now we will compare the projection of these
test samples against the PCs computed for the training samples.

*Now repeat previous steps from [pca.md](pca.md) to load the PCs into
the R environment.*

Here we create the same plot as before, but here we overlay the plot
with the test samples (), showing the expert-provided population
labels. 

```R
pc.test           <- cbind(data.frame(id = rownames(pc.test)),pc.test)
rownames(pc.test) <- NULL
pc.test           <- transform(merge(panel,pc.test),
                               id = as.character(id))
print(plotpca(pc,1,2,dat.more = pc.test))
```

This demonstrates that we can use the PCs to make *predictions* about
any unseen genotype sample as long as (1) we have genotype data for
the same set of SNPs as the training samples, and (2) the genotypes
are encoded in the same way.

As before, try different combinations of PCs.
<commit_msg>Updated first part of Episode 3.<commit_after># Analysis of Genetic Data 1:<br>Making predictions using PCA

Here, we further explore fhe potential of PCA for genetic data by
computing the "projection" of a new set of genotype samples (a "test
set") onto the PCs computed previously from a different set samples
(the "training set").

###Encoding the test genotypes

In order to be able to apply our PCA results to a new set of genotype
samples, we need to make sure that they are represented, or *encoded*,
in the genotype matrix *in the exact same way.* *This can be difficult
to achieve, particularly when the genotype data are obtained using a
different genotyping technology, or in different experimental
conditions.* Specifically, we need to make sure that a "0" genotype in
the training samples means the same thing as a "0" genotype in the
test samples, and likewise for genotypes "1" and "2". In this example,
we can enforce the same encoding this using the `--recode-allele`
option in PLINK:

```bash
module load plink/1.90
cut -f 2,5 1kg_train.bim > alleles.txt
plink --file 1kg_test --recode A-transpose spacex \
  --recode-allele alleles.txt --out 1kg_test
```

We now have the genotype data stored as a matrix with entries that are
0, 1, 2 or missing (NA). Next, we open up R again,

```bash
cd ~/git/genetic-data-analysis-rcc-1/code
R --no-save
```

and in R we enter the following:

```R
# Load the functions for reading data from various files.
source("read.data.R")

# Read the genotypes of the test samples from the .traw file. 
X <- read.geno.traw(fam.file  = "../data/1kg_test.fam",
                    geno.file = "../data/1kg_test.traw")
print(dim(X))

# Read the PCA rotation generated by MATLAB script geno_pca.m.
R <- read.rot.matrix("../results/1kg_train_rot.txt")
print(dim(R))

# Read the mean genotypes outputed by geno_pca.m.
y <- read.mean.genotypes("../results/1kg_train_mean.txt")
```

The genotype matrix we just loaded should have one row for each of the
29 test samples and one column for each of the 156,923 SNPs (same as
the training samples). The rotation matrix is a p x k matrix, where k
is the number of PCs computed.

A simple matrix multiplication completes the projection of the test
samples onto the embedding defined by the PCs.

```R
source("misc.R")
X           <- center.cols(X,y)
X[is.na(X)] <- 0
pc.test     <- X %*% R
```

The only annoying problem was the small number (<1%) of missing
genotypes that make the matrix multiplication less straightforward.
Here we simply set the missing entries to zero after centering the
columns of the matrix. Now we will compare the projection of these
test samples against the PCs computed for the training samples.

*Now repeat previous steps from [pca.md](pca.md) to load the PCs into
the R environment.*

Here we create the same plot as before, but here we overlay the plot
with the test samples (), showing the expert-provided population
labels. 

```R
pc.test           <- cbind(data.frame(id = rownames(pc.test)),pc.test)
rownames(pc.test) <- NULL
pc.test           <- transform(merge(panel,pc.test),
                               id = as.character(id))
print(plotpca(pc,1,2,dat.more = pc.test))
```

This demonstrates that we can use the PCs to make *predictions* about
any unseen genotype sample as long as (1) we have genotype data for
the same set of SNPs as the training samples, and (2) the genotypes
are encoded in the same way.

As before, try different combinations of PCs.

{'fraction_non_alphanumeric': 0.08917865707434053, 'fraction_numerical': 0.010641486810551559, 'mean_word_length': 3.6243936243936243, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4301940', 'n_tokens_mistral': 2091, 'n_tokens_neox': 1915, 'n_words': 952}