# Introduction to BLUPF90 software suite
Docente: Prof. Dra. Daniela Lourenço
Date: 16-17/11/2022
Discente: Larissa Graciano Braga

**BLUPF90 é uma coleção de programas**
Programas para calcular a variância para depois fazer o BLUP (Equação de modelos mistos)

Fortran: melhor linguagem para cálculo -> ainda é utilizado e está sendo "atualizado"

### Main developers
Misztal - criou
Shogo Tsuruta - veio depois do misztal e foi responsável por muitos programas da família
Andres Legarra
Ignacio Aguilar: UY
Yutaka Masuda: atualmente melhor prog
Matias Bermann

## Programas
blupf90 - resolve eq de modelos mistos
remlf90 - resolve covariâncias
airemlf90 - average (libera standar erros, mais usado que EM)
gibbsXf90 - componentes var usando bayesiana
thrgibbsXf90 - para caract não linear, para modelos de threshold

Genomic:
preGS e qc: controle de qualidade, qcf90 é mais eficiente
postGSf - estima valor SNP e faz GWAS
predf90 - pega o efeito de SNP e calcula valor gen genomico
seekparent - checa parentesco, procura pais genotipados para filhos com divergência o parentesco
predictf90 - calcula fenótipo ajustado pelos efeitos fixos

> Input
Todos precisam de arquivo de parâmetros, exceto qcf90
Renumerar os efeitos de 1 até n

- Atualizações (2022?):
blupf90+: junção do blupf90, remlf90, airemlf90
gibbsf90+: junção dos gibbsXf90 e thrgibbsXf90

# Dia 1
## RENUMF90
Renumbers data and pedigree
Creates a par file for blupf90 family
Traces back pedigree for individuals in the data
Para cada animal com fenótipo ou genótipo e rastreia 3 gerações para trás
number of generations (default): 3 back
Performs comprehensive pedigree checks - checa se id de dam e sire é igual
Xref file
Xref id - nova id para os animais genotipados e outra column com
Não renumera arquivo de SNP
Computes inbreeding by default (v>= 1.157)

- Suports
Multiple traits, cria pai fantasma, calc covariâncias

- Input
Arquivo txt separado por espaço, não pode ter comentário (#)
Se tiver info desconhecida sobre dam e sire deve colocar 0 (00 é lido como id de algum animal)
Arquivo parâmetro: características
palavras-chave: DATAFILE (em baixo coloca o nome do arquivo de dados), TRAITS (qual coluna esta a caracteristica), FIELDS_PASSED TO OUTPUT (numero de colunas que vao ser printadas no output), WEIGHT(S) (pesos para var residual - quase n usa isso), .,
EFFECT: uma palavra effect para cada efeito do seu modelo, coloca o numero da coluna que esta o efeito. Ai vc coloca junto: cross (alpha ou number) ou cov
NESTED (aninhar covariavel):coloca a coluna do efeito e também coloca: alpha ou number

- Output
inbreding,
arquivo de dados .dat
arquivo de pedigree renaddxx.ped (renumbered pedigree statistic - xx é o numero do efeito)
arquivo de parâmetros

> exemplo modelo fixo:
DATAFILE
data1.txt
TRAITS
5
FIELDS_PSASSED TO OUTPUT
2 3
WEIGHT(S)
RESIDUAL VARIANCE
1.0
EFFECT
2 cross alpha
EFFECT
3 cross number
EFFECT
4 cov

RANDOM: diagonal (qnd n é correlacionado) ou animal (efeito aleatorio correlacionado entre animais)
OPTIONAL: pe (E permanente) ou mat (maternal) ou mpe (efeito permanente materno - only if mat is used)

keywords afeter RANDOM
FILE: nome do arquivo de pedigree
FILE_POS: ani sire dam alternate_dam yob
SNP_FILE: nome do arquivo de SNP (formato do arquivo de SNP: ID 0111221212122)

PED_DEPTH: default 3 (0 para usar todo mundo, mas nao recomenda usar pelo custo computacional)
GEN_INT: ela n usa, prediz idade dos pais com base do intervalo geração
REC_SEX: fenótipo limitados pelo sexo
UPG_TYPE: yob, in_pedigrees, groups_unisex, group_sex
para pais fantasma (metafounders) - se n tiver pais, vai regredir o VGE para a base da população (inicio do MGA)
Para estimar metafounders: min 1000 fenotipos para cada metafounder - nao usar se banco de dados pequeno
Designar pais com base no ano
INBREEDIN: pedigree, file

Six methods: opção 1 default, opção 6 é a mais rápida (para um numero MUITO grande de dados)

RANDOM_REGRESSION: data
RR_POSITION

(CO)VARIANCES: valor das cov (matriz nxn)
direto | cov direto, materno
cov direto, materno | materno

(CO)VARIANCES
...
COMBINE: combina coluna, deve ser a primeira palavra-chave
7 2 3 4 - coluna nova 7, combinando as col 2 3 e 4

OPTION sol_se (antigo)
OPTION use_yams (recomenda, mais rápido!!!)

Exemplo modelo mistos:
> EFFECT
2 cross alpha
RANDOM
(CO)VARIANCES
0.5
EFFECT # quarto efeito: animal
1 cross alpha
RANDOM
animal
FILE
ped1.txt
FILE_POS
1 2 3 0 0 0
(CO)VARIANCES
0.5

renadd1.ped
colums 4: código com numero grande pra ser usado depois - n preocupe com essa coluna

Problemas: erros de digitação, cov quadrática (crie uma coluna nova com o valor quadratico)
Quando usar renumf90: sempre

Como comparar modelos:
Começa com renumf90 pelo modelos mais complexos
E depois vc retira os efeitos e faz os modelos mais simples
Se quiser retirar animais: retira a linha depois do renumf90

- Quick trick
`renumf90 --help`
`renumf90 --show_template`

## BLUPF90+
combina bulp90, aireml e reml
Default: MME

Se quiser calcular variancia: use `OPTION method VCE`
Gera erro-padrão
Para rodar EM REML - expectation maximization
`OPTION method VCE`
`OPTION EM-REML -1000` (*coloque um numero alto negativo - só para ele soltar valor EM, se colocar positivo ele solta AI*)

Soluções por vários métodos:
PCG: (tenta achar matriz similar a inversa - diagonal da left-hand side)
Para modelo de multicaract use OPTION blksize X (X é o numero de caract do modelo)
SOR:
Inversa pela decomposição de Cholesky: YAMS (reocmendável) ou FSPAK
OPTION use_yams

- arquivo parametro gerado pelo RENUMf90
RANDOM_GROUP
num de efeitos
RANDOM_TYPE
add_an_upginb
Outras opções para dados perdidos (dados, n é pedigree): OPTION missing -999

Calcula reliabilities: `OPTION store_accuracy X` (**X** é o numero do efeito animal)

Nova opção: salva solução com id original no arquivo solutions.original
`OPTION origID`

Problema: NaN tem arquivo pedigree ou data errado

- 2 algoritmos REML: AI e EM

EM REML (iterativo):
1. Coloca valor inicial de variancia,
2. Computa sol
3. Atualiza componentes de variancia (C^uu)
4. Volta para 1

AI-REML (também iterativo)
Estima erro padrão das estimativas e herdabilidades
Alto custo computacional, mais rápido que EM
Compensa começar em com EM e depois passar para A`I OPTION method VCF`; `OPTION EM-REML 10`)

> Problema:
-N converge
Chuta valores muito maiores ou muito menores
-N converge com AI, mas converge com EM
Rode EM com valor menor

blupf90+ --help

## se_covar_function
https://nce.ads.uga.edu/wiki/doku.php?id=readme.aireml
`OPTION se_covar_function label G_eff1_eff2_trt1_trt2 e R_trt1_trt2`
- EXEMPLO (a id do animal esta na col 4)

TRAITS #coluna característica
3 4 5
...
EFFECT
2 cross alpha #EFEITO DE GC # numero da coluna
EFFECT 
4 cross alpha #efeito fixo IPP # tudo eh pela coluna, 
EFFECT
5 cross alpha #efeito fixo PAC
EFFECT
1 cross alpha #animal, o quarto efeito
RANDOM
animal
... 
OPTION se_covar_function H2ipp G_4_4_1_1/(G_4_4_1_1+R_1_1) #herdabilidade posição do animal na 4 (ver efeito do modelo), 1 caracteristica (olhar o TRAITS, se duas: 2_2)
#OPTION se_covar_function RG_ipp_carS G_4_4_1_2/(G_4_4_1_1*G_4_4_2_2)**0.5 #  4 eh o animal (ver efeito modelo), 1_2 correl entre caract 1 e 2 (1 comando para cada tipo de correl X e Y, X e W, Y e W), var dar caract 1 (G_4_4_1_1) vzs var da caract 2 (G_4_4_2_2) (olhar o TRAITS - eu defini que a caract 1 estah na col2, a caract 2 estah na col4 e a caract 3 estah na col5
#OPTION se_covar_function RG_ipp_carG G_4_4_1_3/(G_4_4_1_1*G_4_4_3_3)**0.5

>Examples
OPTION se_covar_function P G_2_2_1_1+G_2_3_1_1+G_3_3_1_1+G_4_4_1_1+R_1_1
OPTION se_covar_function H2d G_2_2_1_1/(G_2_2_1_1+G_2_3_1_1+G_3_3_1_1+G_4_4_1_1+R_1_1)
OPTION se_covar_function H2t (G_2_2_1_1+1.5*G_2_3_1_1+0.5*G_3_3_1_1)/(G_2_2_1_1+G_2_3_1_1+G_3_3_1_1+G_4_4_1_1+R_1_1)
OPTION se_covar_function rg12 G_2_2_1_2/(G_2_2_1_1*G_2_2_2_2)**0.5
-The first function calculates the SD for the total variance for a maternal model with permanent maternal effect, where 2 and 3 are the effect number for the direct and maternal additive genetic effects respectively, and 4 is the effect number for the maternal permanent random effect.
-The second function calculates the heritability for the direct component.
-The third function the total heritability.
-The fourth function calculates the SD of the genetic correlation between traits 1 and 2 for the direct genetic effect (effect number 2)

## gibbsf90+
gibbsXf90
thrgibbsXf90

Numerar categoria de 1 a 5
OPTION cat 0 2 5

`gibbsf90+ parfile.par --samples i --burnin j --interval k`
Burnin: só calcula estatisticas após j amostras, porque a estimativa varia muito (n dá pra confiar)
Interval: a cada 10 amostras, grava 1

gibbsf90: estima covariancias, primeiro vc coloca burnin 0, daí olha o gráfico e decide até quando vai burnin
postgibbsf90: coloca valor burnin
Roda novo gibbsf90+ com novos componentes variância para estimar VG

Modelo Linear ou de treshold? Daniela prefere linear
Daniela: Pode estimar variancias pelo gibbs e rodar o BLUP, n tem problema

Opção para evitar perda por queda de energia: OPTION save_halfway_samples n

Quick trick:
`gibbsf90+ --help`
gibbsf90+ para genomic data:
Rode renumf90 com OPTION animal_order genotypes
Rode gibbsf90+ com OPTION

Problemas comuns:
posição errada, formato
Mispelling, (co)variance matrices not symmetric

output printed on the screen is not saved to any file: you should use direction or pipes to store output
renumf90
rnumf90 renum.paz | tee .log

Background - rodar em bckground

Matriz Z:
Z = M - 2p
p é a MAF; Z é a matriz de genótipo id 012120102102

G = ZZ'/ 2Ep(1-p)

ZZ'= pega a Z e multiplica pela transposta para a gente ter o relacionamento genômico entre animais ("
soma dos SNP")

Relacionamento genomico: cov entre animais
Genomica dá relacionamento observado entre inidividuos
Se LD for alto, vou conseguir uma boa informação com os SNP. Como se usa muitos SNP, temos bons aproximadores do parentesco verdadeiro
A - identico por descendencia
G - identico por estado

Precisamos inverter G
Mas n consegue quando o número de animais > numero de SNP (porque tem dependências lineares)
Como resolver: blending (para aumentar a diagonal) combino a matriz G com uma proporção da matrix A ou com uma identidade
G*= 0.95G + 0.05A
G* = 0.95 + 0.05I
alpha beta
Blending é equivalente ao resíduo poligênico

#### NEM TODOS OS ANIMAIS SÃO GENOTIPADOS
Se fosse, era só pra usar o GBLUP. Mas não é assim!

Podemos dividir a A (matriz de parentesco dos animais)
a11 (n genotipados) a12 (genotipados e nao genotipados)
a21 (igual a a12) a22 (genotipados)
Vc pode substituir/combinar a22 com G -> H

Single step GBLUP (ssGBLUP) usa a H

- Correlação entre G e A
Deve estar entre 0.7-0.9 (mais ou menos do que isso deve ser erro - no pedigree ou nos seus dados)
Importante para saber se incluir genômica vale a pena ou não. Ou para detectar erros

Na H inversa subtrai a a22 pra evitar duplicidade de informação

Tunning: colocar G e a22 na mesma base.
O blupf90: corrige G para ter a mesma média de elementos e diagonal para a22

## SNP QC
preGSf90 faz controle de qualidade genômica (blupf90+ faz tbm, mas n salva arquivo) e cria e inverte A e G
Computa stats
H-1 = A-1 + 0 0
0 G-1a22
Preds salva arquivo "GimAssi" que contem G-1a22

`OPTION SNP_FILE marker.geno` (no renumf90 SNP_FILE e em baixo marker.geno)

`OPTION map_file file.map`
Formato: com header

O predsf90 faz QC e inverte as matrizes G e A automaticamente (demanda tempo e processamento). Se vc quiser só fazer QC use algumas options

OPTION whichfreq x
1 0.5 (p/ metafounders)

- Tuning
OPTION tunedG 2

- Population structure
OPTION plot_pca
OPTION extra_info_pca
(para raças diferentes)

Constrói a G usando VanRaden (2008), a22 usando Colleau (2002) e A-1 pelo método henderson e quatis(?)

Atividade
`ulimit -s unlimited` # para liberar memória

curl http://nce.ads.uga.edu/wiki/lib/exe/fetch.php?media=lab1_ufv.zip -o lab1.zip
renumf90 --show-template

#montar cartão de parâmetros usando os arquivos do lab1 (nao colocar espaços, excepcionalmente apos FIELDS.. e WEIGHT(S))

DATAFILE
data3.txt
TRAITS
4
FIELDS_PASSED TO OUTPUT
1
WEIGHT(S)

RESIDUAL_VARIANCE
0.6
EFFECT
3 cross alpha # efeito de sexo
EFFECT
1 cross alpha # efeito de animal
RANDOM  # estou informando que a col 1 é o efeito aleatório 
animal
FILE
ped3.txt
FILE_POS
1 2 3 0 0
SNP_FILE
snp3.2k
PED_DEPTH
3
(CO)VARIANCES
0.4 # dado no exercicio, na vrdd e so a variância pq é só uma caract

Linha de comando: `renumf90 my.par | tee renumf90.log &`

- Exercicio 3
What is the content of snp3.2k_XrefID?
Nova id e id original dos animais que foram genotipados

- Exercício 4
Uso o renf.par para rodar o preGSf90

# BLUPF90 parameter file created by RENUMF90
DATAFILE
 renf90.dat
NUMBER_OF_TRAITS
           1
NUMBER_OF_EFFECTS
           2
OBSERVATION(S)
    1
WEIGHT(S)

EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT[EFFECT NESTED]
 2         2 cross
 3     12010 cross
RANDOM_RESIDUAL VALUES
  0.60000
 RANDOM_GROUP
     2
 RANDOM_TYPE
 add_an_upginb
 FILE
renadd02.ped
(CO)VARIANCES
  0.40000
OPTION SNP_file snp3.2k
OPTION map_file mrkmap.txt
OPTION saveCleanSNPs
OPTION msg 100

Linha de comando: `preGSf90 renf90.par | tee preGSf90.log &`

Which quality checks for both SNP and animals were
done by default?
Call rate, MAF, Monomorphic SNP, mendelian conflicts
> Quality Control - SNPs with Call Rate < callrate ( 0.90) will removed: 0
Quality Control - SNPs with MAF < minfreq ( 0.05) will removed: 82
Quality Control - Monomorphic SNPs will be removed: 0
Quality Control - Removed Animals with Call rate < callrate ( 0.90): 0
Quality Control - Check Parent-Progeny Mendelian conflicts

Are there any duplicated genotypes?
2 animais: 2024 (old ID = 2966) e 2020 (old ID =8818)
> Possible genotype duplicates samples
i-j number of sample , i-j renumber Id, G(i,j), G(i,i), G(j,j), r(i,j)
2024 2020 2966 8818 1.0653 1.0653 1.0653 1.0000

What is the correlation between G and A22?
> Correlation all elements G & A 0.780

Check averages of G and A22
G: 0.992 e A22: 1.02

> Statistic of Genomic Matrix
N Mean Min Max Var
Diagonal 2024 0.992 0.816 1.358 0.004
Off-diagonal 4094552 -0.000 -0.234 1.065 0.007

- Exercicio 5

# BLUPF90 parameter file created by RENUMF90
DATAFILE
 renf90.dat
NUMBER_OF_TRAITS
           1
NUMBER_OF_EFFECTS
           2
OBSERVATION(S)
    1
WEIGHT(S)

EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT[EFFECT NESTED]
 2         2 cross
 3     12010 cross
RANDOM_RESIDUAL VALUES
  0.60901
 RANDOM_GROUP
     2
 RANDOM_TYPE
 add_an_upginb
 FILE
renadd02.ped
(CO)VARIANCES
  0.36197
#OPTION map_file mrkmap.txt
#OPTION SNP_file snp3.2k
#OPTION saveCleanSNPs
#OPTION no_quality_control
OPTION createGInverse 0
OPTION createA22Inverse 0
OPTION createGimA22i 0
OPTION use_yams
OPTION method VCE
OPTION se_covar_function H2_1 G_2_2_1_1/(G_2_2_1_1+R_1_1)
OPTION sol se
#OPTION EM-REML -200 # para rodar EM-REML retirar hashtag, o valor deve ser negativo mesmo
OPTION origID
OPTION msg 100
#OPTION store_accuracy 2

Estimativa de tempo que rodou `time blupf90+ renum.par | tee blupf90_BLUP_AIREML.log`

Tempo para BLUP AIREML

 * FINISHED (BLUPF90): 11-17-2022  13h 18m 41s 612

real    0m2.253s
user    0m11.264s
sys     0m1.071s

Verificar as variâncias e rodar novamente com os valores corrigidos

Final Estimates
 Genetic variance(s) for effect  2
  0.36197
 Residual variance(s)
  0.60901

Após, corrigir os campos de `RANDOM_RESIDUAL VALUES` com 0.60901 e `(CO)VARIANCES` com 0.36197

Rodar para EM-REML
Como? retirando hashtag da option em-reml, verificar o número de rodadas, se você colocou 20 e o programas iterou 22 vezes, você vai precisar aumentar o valor de iterações - no exemplo da aula, primeiro colocaram 20 ai o blup rodou 26 vezes, dai foram ajustando o numero de iterações para EM-REML até ficar suficiente

Como vejo o round?

 -2logL =    26900.7548501640       : AIC =    26904.7548501640
  In round            1  convergence=  9.069627257756786E-013
 round            1
 new r
  0.60901
 new G
  0.36197

Eu quero ver - situação ideal: `At round: 1 Converge in fewer rounds than EM-REML rounds: 200`

Rodando o EM-REML: `time blupf90+ renf90.par | tee blupf90_BLUP_EM-REML.log`

 * FINISHED (BLUPF90): 11-17-2022  14h 00m 36s 518

real    0m2.308s

- Exercicio 6
Retirar a hashtag da OPTION store_accuracy 2
Fazer isso para ssGBLUP (retirar hashtag relacionadas a SNP), e BLUP (options relacionadas a snp com #)

ssGBLUP

# BLUPF90 parameter file created by RENUMF90
DATAFILE
 renf90.dat
NUMBER_OF_TRAITS
           1
NUMBER_OF_EFFECTS
           2
OBSERVATION(S)
    1
WEIGHT(S)

EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT[EFFECT NESTED]
 2         2 cross
 3     12010 cross
RANDOM_RESIDUAL VALUES
  0.60000 
 RANDOM_GROUP
     2
 RANDOM_TYPE
 add_an_upginb
 FILE
renadd02.ped
(CO)VARIANCES
  0.40000
OPTION SNP_file snp3.2k_clean
OPTION no_quality_control
#OPTION saveCleanSNPs
OPTION use_yams
#OPTION method VCE
OPTION se_covar_function H2_1 G_2_2_1_1/(G_2_2_1_1+R_1_1)
#OPTION EM-REML -200
OPTION origID
OPTION msg 100
OPTION sol se
OPTION store_accuracy 2

Acurácias no arquivo `acc_bf90`

Como ler o solutions:

trait effect level(new id) original solution s.e
1    1    1    media_para_efeito1    se_mediano_para_efeito1
1    2    2    media_para_efeito1    se_mediano_para_efeito1
1    1    1    animail 1              se para o animal 1
...  n    n    animal n               se para o animal n

# Dia 2
## Validação
Para comparar modelos! e ver quão benéfica a seleção genômica é
Método estatíco: validação por Prediction (validation) accuracy
Aqui temos um valor para toda a população
Sinonimos: vários

Theoretical accuracy - 1 acurácia de cada animal
Sinonimo: reliability

### Theoretical accuracy
Fala quanto os dados vão mudar após incluir mais info. É uma por indivíduo e é baseada no modelo

Baseada no PEV: variância do erro de predição (na formula do pev tem standard erro)

Aumenta quando a gnt aumenta a quantdade de info de um animal

Acurácia: 0.05. Possible change: 2.49
Possible change: intervalo de confiança 95% = `EBV +- 1.96 * standard error`
Reliability != BIF accuracy != accuracy -> muda raiz quadrada e 1- na fórmula delas
BIF é a que penaliza mais
Leite: a maioria dos touros tem muitas filhas, aí para que muitos animais não tenham alta acurácia, eles não usam a raiz quadrada na fórmula

### Prediction accuracy
Correlação entre TBV e EBV.
Um valor único para a população
A gnt calcula acurácia com e sem genômica pra ver se vale a pena

Accurary = cor(benchmark, u^)
Benchmarl: ponto de referência que pode ser progeny yield deviation (leite), deregressed EBV (leite, boars), progeny yield deviation (boars), adajusted phenotypes (beef, pig, aves ...)

Como: grupo training e grupo validation

Benchmark: uso toda info de pedigree, P e G
Para calc GEBV e EBV:reduced data, training data (inicio-intervalo de geração) e validation data (intervalo de geração-ultimo dado)
ex: inicio-2013 e 2013-2017
Macho tem que ter, pelo menos, 10 filhas nos dados completos e, no mín, 100 touros

dEBVcomplete = b0 + b1 * EBVreduzido
dEBVcomplete = b0 + b1 * GEBVreduzido
b0= bias, viés de nível; b1: viés de dispersão; R2
nós esperamos que: R2(GBEV) > R2(EBV); b0(GBEV) e b0(EBV) qual é mais próximo de 0; b1(GBEV) e b1(EBV) qual é mais próximo de 1
Se b1<1, meu (G)EBV está inflacionado. Se b1>1, meu (G)EBV está deflacionado

Dispersão: olhar dispersão dos (G)EBV dos animais jovens pq se estiverem muito dispersos ou pouco dispersos a sua "régua" de seleção pode incluir ou excluir animais, respectivamente, de forma errada

Adjusted phenotypes benchmark
Predictive ability of (g)EBV = cor (Yadj, (G)EBVreduced)

- Linear Regression metrics
paper: Legarra e Reverter
New method: está sendo mais usado
Reduced and complete data
Mesma divisão dos dados que os métodos anteriores (validation animals have phenotypes in the complete data but not in the reduced data)
Benchmark: complete (G)EBV
Compares EBV with EBV e GEBV com GEBV
Como? vai lá no renf90 e corta os dados de 4 anos atrás e procura os touros que podem participar dessa validação
A acurácia que eles estão usando é a mais simples

Se Cor(GEBVcompleto, GEBVreduzido) = 0.80 e Cor(EBVcompleto, EBVreduzido) = 0.58 -> o GEBV é um bom preditor

- Exercício 7
GEBV:
cd ssGBLUP
Linha de comando: `renumf90 my_ssGBLUP.par | tee renumf90_ssGBLUP.log &`

DATAFILE
../data3.txt
TRAITS
4
FIELDS_PASSED TO OUTPUT
1 2
WEIGHT(S)

RESIDUAL_VARIANCE
0.6
EFFECT
3 cross alpha
EFFECT
1 cross alpha
RANDOM
animal
FILE
../ped3.txt
FILE_POS
1 2 3 0 0
SNP_FILE
../snp3.2k
PED_DEPTH
3
(CO)VARIANCES
0.4

`awk '$4!=5' ../renf90.dat > renf90.dat.reduced` # coloquei 5 pq tem uma coluna a mais no meu (a de id original)

Cartão reduzido

# BLUPF90 parameter file created by RENUMF90
DATAFILE
 renf90.dat.reduced
NUMBER_OF_TRAITS
           1
NUMBER_OF_EFFECTS
           2
OBSERVATION(S)
    1
WEIGHT(S)

EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT[EFFECT NESTED]
 2         2 cross    # efeito de sexo
 3     12010 cross    # efeito de animal
RANDOM_RESIDUAL VALUES
  0.61222       # valor calculado nos exercicios anteriores
 RANDOM_GROUP
     2           # o efeito aleatorio é o segundo efeito
 RANDOM_TYPE
 add_an_upginb
 FILE
renadd02.ped
(CO)VARIANCES
  0.35532       # valor calculado nos exercicios anteriores
OPTION SNP_file ../snp3.2k
OPTION use_yams
#OPTION method VCE
OPTION se_covar_function H2_1 G_2_2_1_1/(G_2_2_1_1+R_1_1)
#OPTION EM-REML -200
OPTION origID
OPTION msg 100
OPTION sol se
OPTION store_accuracy 2 # vai guardar o efeito 2

`blupf90+ renf90.par | tee blupf90_ssGBLUP_partial.log`

Mudando o nome do arquivo para manter separado (n sobrescrever)
`mv solutions_ssGLUP_reduced.orig solutions_ssGLUP_partial.orig`

- EBV
Cartão de paramentros para BLUPF90 para EBV, sem genômica

# BLUPF90 parameter file created by RENUMF90
DATAFILE
 renf90.dat.reduced
NUMBER_OF_TRAITS
           1
NUMBER_OF_EFFECTS
           2
OBSERVATION(S)
    1
WEIGHT(S)

EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT[EFFECT NESTED]
 2         2 cross
 3     12010 cross
RANDOM_RESIDUAL VALUES
  0.68990
 RANDOM_GROUP
     2
 RANDOM_TYPE
 add_an_upginb
 FILE
renadd02.ped
(CO)VARIANCES
  0.36200
#OPTION SNP_file ../snp3.2k
#OPTION no_quality_control
OPTION use_yams
#OPTION method VCE
OPTION se_covar_function H2_1 G_2_2_1_1/(G_2_2_1_1+R_1_1)
#OPTION EM-REML -200
OPTION origID
OPTION msg 100
OPTION sol se
OPTION store_accuracy 2

`blupf90+ renf90.par | tee blupf90_BLUP_partial.log`

Mudando o nome para saber
`mv solutions.orig solutions_BLUP_partial.orig`

NÃO FIZ: Calculando a correlação entre os EBV e GEBV dos animais da 5 geração com os seus respectivos TBV

> getwd()

> partial=read.table("partial/solutions.orig", header=T)
> whole=read.table("complete/solutions.orig", header=T)

# Correlação entre partial e completo
> cor.test(partial$solution, whole$solution, method="pearson")

        Pearson's product-moment correlation

data:  partial$solution and whole$solution
t = 315.35, df = 12010, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.9426269 0.9464820
sample estimates:
     cor
0.944587

# Correlação entre TBV e completo apenas para a quinta geração

- Exercício 8
Fazer como na vida real, a gnt nao sabe o TBV

renumf90 roda 1x # só uma vez para comparar modelos, se vc rodar o renum em várias pastas vai dar ruim pq o renum renumera de forma diferente
blupf90 roda 4x
Comparar no LR:
BLUP: EBVpartial, EBVwhole
ssGBLUP: GEBVpartial, GEBVwhole

~~Curiosidade: Cartão de parâmetros do mateus BLUP~~

# BLUPF90 parameter file created by RENUMF90
DATAFILE
 renf90.dat
NUMBER_OF_TRAITS
           1
NUMBER_OF_EFFECTS
           2
OBSERVATION(S)
    1
WEIGHT(S)

EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT[EFFECT NESTED]
 2         2 cross
 3     12010 cross
RANDOM_RESIDUAL VALUES
  0.60000
 RANDOM_GROUP
     2
 RANDOM_TYPE
 add_an_upginb
 FILE
renadd02.ped
(CO)VARIANCES
  0.40000
OPTION SNP_file snp3.2k
OPTION map_file mrkmap.txt
#OPTION saveCleanSNPs
OPTION use_yams
OPTION method VCE
OPTION se_covar_function H2_1 G_2_2_1_1/(G_2_2_1_1+R_1_1)
OPTION EM-REML -200
OPTION origID
OPTION msg 100

## GWAS

O que é a matriz W?

#### SNP-BLUP
Z* a = u
> Matrix Z: animal por SNP
a: efeito de SNP (calculado pelo postGSf90)
u = vetor de DGV
Se Z = 2000 animais * 50k SNPs e a = 50000 * 1
u = 2000 * 1, ou seja, o valor genômico para cada animal

GBLUP e SNP-BLUP são equivalentes pq a variância dos modelos é igual

Slide: SNP effects in ssGBLUP
alfa=blending (0.95 * 0.05, da aula do dia 1)
b = coef de tunning
Regressão para obter valor de SNP

Posso fazer indirect prediction e GWAS com isso aqui de cima
GWAS

Onde colocar a barra de treshold (acima da barra é significante)? Correção de Bonferroni é 0.05/nº SNP
Pq -logPvalue? O valor de log é pq os valores ficam muito pequenos (0.0000002). Também pode colocar a variância explicada pelo SNP

##### Indirect prediction (IP)
DEP-interina

Se for comparar IP e GEBV - precisa colocar na escala de pedigree, para ter a mesma base genética dos animais da avaliação

Rodar blupf90 e salvar os arquivos limpos - OPTION saveGInverse + OPTION
Rodar postGSf90 : option use inverse... e com o mapa de SNP

Output:
snp_sol col 7, peso dos SNP

predf90: arquivo de SNP dos novos animais para fazer a predição, é o mesmo que o arquivo limpo
roda só linha de comando `predf90 --snpfile newgen.txt --use_mu_hat` # mu para ficar na mesma base do GEBV

##### GWAS
Remove estrutura para parentesco, uma regressão para cada SNP

Métodos current (p/ humano) n funciona por causa das nossas limitações
dEBV: projeta P das filhas no touro

ssGWAS!
d = se acima de 1 quer dizer que distancia da normal

Se SNP é verdadeiro, a gnt ve uma trilha de SNP por causa de LD

P-value in ssGWAS
1-2: blupf90 (solta arquivo binário XX_ija que tem os coeficientes para animais genotipados)
3-6: postGS

blupf90
OPTION SNP_FILE clean, MAP_FILE clean, saveGInverse, saveA22Inverse, snp_p_value
postGSf90
option windows_variance X # 1Mb de janela= 50k 20 SNP
Output: arquivo que começa com P são os que tem p-valor, arquivo que começa cm s, snp_sol

GWAS não é causal, serve para aprender um pouco mais sobre a arquitetura genômica

###### Lab 2
- Exercício 1
`renumf90 renum.par | tee renumf90.log`

renf90.par

# BLUPF90 parameter file created by RENUMF90
DATAFILE
 renf90.dat
NUMBER_OF_TRAITS
           1
NUMBER_OF_EFFECTS
           2
OBSERVATION(S)
    1
WEIGHT(S)

EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT[EFFECT NESTED]
 2         2 cross
 3     12010 cross
RANDOM_RESIDUAL VALUES
  0.60000
 RANDOM_GROUP
     2
 RANDOM_TYPE
 add_an_upginb
 FILE
renadd02.ped
(CO)VARIANCES
  0.40000
OPTION SNP_file snp3.2k
OPTION map_file mrkmap.txt
OPTION no_quality_control
OPTION saveGInverse # salva a matriz G inversa
OPTION saveA22Inverse # salva a matriz A22
OPTION snp_p_value # para criar arquivo xx_ija que calcula do p-value - o postGS nao calcula pvalue
OPTION origID 

Rodar BLUP90 com as opções que preciso
`blupf90+ renf90.par | tee blupf90_to_GWAS.log`

Rodar o postGSf90, o cartão de parâmetros:

# BLUPF90 parameter file created by RENUMF90
DATAFILE
 renf90.dat
NUMBER_OF_TRAITS
           1
NUMBER_OF_EFFECTS
           2
OBSERVATION(S)
    1
WEIGHT(S)

EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT[EFFECT NESTED]
 2         2 cross
 3     12010 cross
RANDOM_RESIDUAL VALUES
  0.60000
 RANDOM_GROUP
     2
 RANDOM_TYPE
 add_an_upginb
 FILE
renadd02.ped
(CO)VARIANCES
  0.40000
OPTION SNP_file snp3.2k
OPTION map_file mrkmap.txt
OPTION no_quality_control
OPTION readGInverse
OPTION readA22Inverse
OPTION snp_p_value
OPTION windows_variance 1 # para verificar o pvalor de todos os SNP (n usa windows)

`Rodando: postGSf90 postGSf90.par | tee postGSf90.log`

Gerar o gráfico: `Rscript *.R` no R

Mas eu editei o arquivo, abri o script no R para colocar a treshold line:

a=0.05/nrow(yyy) #formula para calculo do treshold
a=-log10(a) # colocar na escala do grafico

E antes de salvar o gráfico (antes dev.off) coloquei abline

Arquivos para olhar
snp_pred #predição do valor genômico de cada SNP
windows_segments # me dá as genomic windows
windows_variance

Agora eu vou fazer a predição para os animais novos (que não tem fenótipos)
Rodar: `predf90`
Ele pergunta a name of genotype file? eu vou colocar `new_animals`
Output: SNP_predictions # me dá o GEBV nos animais jovens (sem P)