# Filogenômica de *Utricularia* (seção *Utricularia*: Lentibulariaceae): uma possível história de evolução reticulada e estado de conservação de um grupo de plantas carnívoras
### FAPESP #20/11053-1
## Análise de dados a partir de Genotype-by-Sequencing (GBS)
Para a análise dos dados foi utilizado o pipeline ***ipyrad***, sendo que foram feitos testes a partir de montagem *de novo* e com referência. Ademais foram feitos testes com diferentes quantidades de *missing data*, modificando principalmente os parâmetros de quantidades mínimas de *loci* por população. Abaixo encontram-se as comparações mais interessantes para discussão do trabalho:
O comando básico para todas as análises fora realizado a partir do comando:
```python!=
ipyrad -f -p params-Utricularia.txt -s 1234567
```
## Criar "ramificação" da análise somente com as amostras de Utricularia foliosa (27 amostras)
```python=
ipyrad -p params-Utricularia.txt -b Ufoliosa_T1 3336_2 3336_5 3376_1 3376_2 3376_3 3376_4 3376_6 3380_1 3380_2 3380_3 3380_5 pop1_10 pop1_11 pop1_12 pop1_13 pop1_14 pop1_15 pop1_2 pop1_6 pop1_7 pop1_8 pop1_9 pop5_1 pop5_2 pop5_3 pop5_4 pop5_5
```
## Criar "ramificação" da análise somente com as amostras de Utricularia gibba (56 amostras)
```python=
ipyrad -p params-Utricularia.txt -b Ugibba 3301_1 3301_2 3301_3 3301_4 3301_5 3304_2 3307_1 3307_2 3309_2 3309_3 3309_4 3309_5 3309_6 3321_1 3321_2 3321_3 3321_4 3321_5 3322_1 3322_2 3322_3 3322_4 3322_5 3322_6 3322_7 3331_1 3331_2 3331_3 3334_3 3334_4 3335_1 3335_2 3335_3 3335_4 3337_1 3337_2 3337_4 3337_6 3337_7 3367_2 3367_3 3367_4 3367_5 3368_0 3379_1 3379_2 3379_5 3379_6 3382_1 3646_1 3646_2 3646_3 3646_4 3646_5 3646_6 3379_3
```
## Criar "ramificação" da análise somente com as amostras de Utricularia poconensis (5 amostras)
```python=
ipyrad -p params-Utricularia.txt -b Upoconensis 2110_3 2110_1 3349_4 3349_5 3349_2
```
## Análises realizadas em *Utricularia foliosa*
#### Arquivos "fundamentais" e comandos utilizados nas análises
Para utilizar o pipeline do *ipyrad* para este trabalho utilizamos os arquivos:
1. com as referências dos barcodes [barcodes_path](https://drive.google.com/file/d/1IRSwgJgbU35-Y-7ngWW4YJLIwIzUhQfN/view?usp=sharing):
2. com identificação das amostras de acordo com as populações [pop_assign_file](https://drive.google.com/file/d/1CyPJzu1zw_3XlrPHqau0p6f3nDUmRRl5/view?usp=sharing):
```
ipyrad -f -p params-Ufoliosa_T1.txt -s 34567
```
Sendo que as variações dos parâmetros para cada Teste foram feitas no arquivo "params-Ufoliosa_T*N*.txt".
:::info
Haja vista que *U. foliosa* é possivelmente uma espécie poliplóide modificamos o parâmetro [max_alleles_consens] = 2 para que a análise assuma os valores de heterozigozidade calculados pelo pipeline.
[Dados de simulações de trabalhos anteriores ainda não publicados do laboratório de sistemática vegetal a partir de *likelihood* dos números cromossômicos dos ancestrais de *Utricularia* nos permitem inferir que *U.foliosa* deve ser poliplóide (n=21)](https://drive.google.com/file/d/1lw_QrxxQZHKwkJOMcd1PjayRZkK33G2V/view?usp=sharing)
:::
| População | #indivíduos | Código | Localidade |
| -------- | -------- | -------- | -------- |
| 3336 | 2 | MC_SP_APA | APA Várzea, Mogi das Cruzes, SP |
| 3376 | 5 | CG_MG |Campos Gerais, Minas Gerais |
| 3380 | 4 | MC_SP_CR | Cerâmica, Mogi das Cruzes, SP |
| pop1 | 11 | UB_SP | Ubarana, São Paulo |
| pop5/2209 | 5 | MC_SP_BOTU| Botujuru, Mogi das Cruzes, SP |
#### Teste 1: Montagem *de novo* sem limitação de quantidade mínima de *loci*
Para este teste realizamos as análises sem referência, ou seja, com montagem *de novo* e clusterização das *reads*
###### Arquivo params-Ufoliosa_T1.txt:
```python!=
------- ipyrad params file (v.0.9.84)---------------------------------------
Ufoliosa ## [0] [assembly_name]: Assembly name. Used to name output directories for assembly steps
/srv/projects/GBS/Analyses_new ## [1] [project_dir]: Project dir (made in curdir if not present)
/srv/projects/GBS/RAW/EM89_GBS_Saura_S1_L001_R1_001.fastq ## [2] [raw_fastq_path]: Location of raw non-demultiplexed fastq files
/srv/projects/GBS/Barcodes/Barcode_file.csv ## [3] [barcodes_path]: Location of barcodes file
## [4] [sorted_fastq_path]: Location of demultiplexed/sorted fastq files
denovo ## [5] [assembly_method]: Assembly method (denovo, reference)
## [6] [reference_sequence]: Location of reference sequence file
gbs ## [7] [datatype]: Datatype (see docs): rad, gbs, ddrad, etc.
TGCAG, ## [8] [restriction_overhang]: Restriction overhang (cut1,) or (cut1, cut2)
5 ## [9] [max_low_qual_bases]: Max low quality base calls (Q<20) in a read
33 ## [10] [phred_Qscore_offset]: phred Q score offset (33 is default and very standard)
6 ## [11] [mindepth_statistical]: Min depth for statistical base calling
6 ## [12] [mindepth_majrule]: Min depth for majority-rule base calling
10000 ## [13] [maxdepth]: Max cluster depth within samples
0.85 ## [14] [clust_threshold]: Clustering threshold for de novo assembly
0 ## [15] [max_barcode_mismatch]: Max number of allowable mismatches in barcodes
2 ## [16] [filter_adapters]: Filter for adapters/primers (1 or 2=stricter)
35 ## [17] [filter_min_trim_len]: Min length of reads after adapter trim
4 ## [18] [max_alleles_consens]: Max alleles per site in consensus sequences
0.05 ## [19] [max_Ns_consens]: Max N's (uncalled bases) in consensus
0.05 ## [20] [max_Hs_consens]: Max Hs (heterozygotes) in consensus
1 ## [21] [min_samples_locus]: Min # samples per locus for output
0.2 ## [22] [max_SNPs_locus]: Max # SNPs per locus
8 ## [23] [max_Indels_locus]: Max # of indels per locus
0.5 ## [24] [max_shared_Hs_locus]: Max # heterozygous sites per locus
0, 0, 0, 0 ## [25] [trim_reads]: Trim raw read edges (R1>, <R1, R2>, <R2) (see docs)
0, 0, 0, 0 ## [26] [trim_loci]: Trim locus edges (see docs) (R1>, <R1, R2>, <R2)
* ## [27] [output_formats]: Output formats (see docs)
/home/srsilva/GBS/Populations_ref_new/Pop_fileUfoliosa_T1.csv ## [28] [pop_assign_file]: Path to population assignment file
## [29] [reference_as_filter]: Reads mapped to this reference are removed in step 3
```
#### Teste 2: Montagem *de novo* com limitação de quantidade mínima de *loci*
Para este teste limitamos a quantidade de *loci* a serem considerados por população em 4 *loci* e 2 *loci* para a população identificada como <i>3336</i>, já que conseguimos amostrar somente 2 indivíduos dessa população. Com a modificação deste parâmetro restringimos mais a análise, permitindo a construção de matrizes com menos dados faltantes e possívelmente menor quantidade de *loci* com potencial de serem falsos positivos, já que ele deverá aparecer em pelo menos 4 indivíduos por população.
###### Arquivo params-Ufoliosa_T2.txt:
Para agregar o parâmetro acima, fizemos a seguinte modificação na última linha do arquivo [pop_assign_file]:
```python!=
# pop1:2 pop2:4 pop3:4 pop4:4 pop5:4
```
#### Outros testes
Também fora feito o teste sem a população 3336, haja vista que só há 2 indivíduos sequenciados até o momento, sontudo a retirada dessa população não permitiu grande aumento do resgate de loci das demais amostras, portanto optou-se por manter essa população para as análises atuais, sendo que serão realizados mais extrações de indivíduos desta mesma população para incorporar na análise numa próxima etapa.
#### Teste 3: Montagem **com referência** sem limitação de quantidade mínima de *loci*
:::danger
Para as análises com genoma referência, utilizamos a montagem dos genomas nuclear (artigo em desenvolvimento) e cloropastidial (publicado em Silva et al., 2017) que foi feita a partir de sequenciamento *genome skimming* de dados produzidos no projeto FAPESP #18/02285-6
:::
###### Arquivo params-Ufoliosa_T3.txt:
Para agregar o genoma referência, modificamos os parâmetros no arquivo **param-Ufoliosa_T3.txt**:
```python!=
reference ## [5] [assembly_method]: Assembly method (denovo, reference)
/srv/projects/GBS/reference/Ufoliosa/Ufoliosa_nuclear.fas ## [6] [reference_sequence]: Location of reference sequence file
```
#### Teste 4: Montagem **com referência** com limitação de quantidade mínima de *loci*
Para agregar o parâmetro acima, fizemos a seguinte modificação na última linha do arquivo [pop_assign_file]:
```python!=
# pop1:2 pop2:4 pop3:4 pop4:4 pop5:4
```
###### Arquivo params-Ufoliosa_T4.txt:
##### Resultados da análise
## Análises realizadas em *Utricularia gibba*
Para as análises com referências, foram feitos testes tanto com o genoma nuclear de Ibarra-Laclette et al. (2013), que foi montado com *reads* curtos, como o genoma de Lan et al. (2017), que fora montado com *reads* longos. Os melhores resultados foram obtidos com o genoma de Lan et al. (2017), sendo estes os apresentados abaixo:
| População | Número de indivíduos| Código | Localidade |
| -------- | -------- | -------- | -------- |
| 3301| 5|AT/MS|Aparecida do Taboado (MS)|
| 3304| 1|CH/GO|Chapadão do Céu (GO)|
| 3307| 2|POC_BR/MT|Poconé (MT) - Barrinha|
| 3309| 5|SON/MS|Sonora (MS)|
| 3321| 5|POC_PI/MT|Poconé (MT) - Pixaim|
| 3322| 6|AQUI/MS|Aquidauana (MS)|
| 3331| 3|NA/MS|Nova Andradina (MS)|
| 3334| 2|PL/SP|Planalto (SP)|
| 3335| 4|MC_APA/SP|Mogi das Cruzes (SP) - APA Várzea do Rio Tietê|
| 3337| 5|SA/SP|Sales (SP)|
| 3367| 4|CA/MG|Carrancas (MG) - Cachoeira da Zilda|
| 3368| 1|CA2/MG|Carrancas (MG) - Cascata do Índio|
| 3379| 5|MC_CR/SP|Mogi das Cruzes (SP) - Cerâmica|
| 3382| 1|BI/SP|Biritiba Mirim (SP), entrada da trilha para a Pedra do Garrafão|
| 3646| 6|SA2/SP|Sales (SP)|
#### Arquivos "fundamentais" e comandos utilizados nas análises
Para utilizar o pipeline do *ipyrad* para este trabalho utilizamos os arquivos:
1. com as referências dos barcodes [barcodes_path](https://drive.google.com/file/d/1IRSwgJgbU35-Y-7ngWW4YJLIwIzUhQfN/view?usp=sharing):
2. com identificação das amostras de acordo com as populações [pop_assign_file](https://drive.google.com/file/d/1abWfCRYVvplrnIv_h14U6oT1qjUzd9e3/view?usp=sharing):
O comando básico para todas as análises fora realizado a partir do comando:
```python!=
ipyrad -f -p params-Ugibba.txt -s 1234567
```
Sendo que as variações dos parâmetros para cada Teste foram feitas no arquivo "params-Ugibba_T*N*.txt".
:::info
Haja vista que *U. gibba* é possivelmente uma espécie diplóide ([Ibarra-Laclette et al., 2013](https://www.nature.com/articles/nature12132)) modificamos o parâmetro [max_alleles_consens] = 2 para que a análise assuma os valores de heterozigozidade calculados pelo pipeline.
:::
#### Teste 1: Montagem *de novo* sem limitação de quantidade mínima de *loci*
Para este teste realizamos as análises sem referência, ou seja, com montagem *de novo* e clusterização das *reads*
###### Arquivo params-Ugibba_T1.txt:
```python!=
------- ipyrad params file (v.0.9.84)-------------------------------------------
Ugibba ## [0] [assembly_name]: Assembly name. Used to name output directories for assembly steps
/srv/projects/GBS/Analyses ## [1] [project_dir]: Project dir (made in curdir if not present)
/srv/projects/GBS/RAW/EM89_GBS_Saura_S1_L001_R1_001.fastq ## [2] [raw_fastq_path]: Location of raw non-demultiplexed fastq files
/srv/projects/GBS/Barcodes/Barcode_file.csv ## [3] [barcodes_path]: Location of barcodes file
## [4] [sorted_fastq_path]: Location of demultiplexed/sorted fastq files
denovo ## [5] [assembly_method]: Assembly method (denovo, reference)
## [6] [reference_sequence]: Location of reference sequence file
gbs ## [7] [datatype]: Datatype (see docs): rad, gbs, ddrad, etc.
TGCAG, ## [8] [restriction_overhang]: Restriction overhang (cut1,) or (cut1, cut2)
5 ## [9] [max_low_qual_bases]: Max low quality base calls (Q<20) in a read
33 ## [10] [phred_Qscore_offset]: phred Q score offset (33 is default and very standard)
6 ## [11] [mindepth_statistical]: Min depth for statistical base calling
6 ## [12] [mindepth_majrule]: Min depth for majority-rule base calling
10000 ## [13] [maxdepth]: Max cluster depth within samples
0.85 ## [14] [clust_threshold]: Clustering threshold for de novo assembly
0 ## [15] [max_barcode_mismatch]: Max number of allowable mismatches in barcodes
2 ## [16] [filter_adapters]: Filter for adapters/primers (1 or 2=stricter)
35 ## [17] [filter_min_trim_len]: Min length of reads after adapter trim
2 ## [18] [max_alleles_consens]: Max alleles per site in consensus sequences
0.05 ## [19] [max_Ns_consens]: Max N's (uncalled bases) in consensus
0.05 ## [20] [max_Hs_consens]: Max Hs (heterozygotes) in consensus
1 ## [21] [min_samples_locus]: Min # samples per locus for output
0.2 ## [22] [max_SNPs_locus]: Max # SNPs per locus
8 ## [23] [max_Indels_locus]: Max # of indels per locus
0.5 ## [24] [max_shared_Hs_locus]: Max # heterozygous sites per locus
0, 0, 0, 0 ## [25] [trim_reads]: Trim raw read edges (R1>, <R1, R2>, <R2) (see docs)
0, 0, 0, 0 ## [26] [trim_loci]: Trim locus edges (see docs) (R1>, <R1, R2>, <R2)
* ## [27] [output_formats]: Output formats (see docs)
## /home/srsilva/GBS/Populations_ref/Pop_fileUgibba.csv ## [28] [pop_assign_file]: Path to population assignment file
## [29] [reference_as_filter]: Reads mapped to this reference are removed in step 3
```
###### Resultados de Ugibba_T1:
O programa ipyrad gera um relatório resumido dos resultados [Ugibba_T1_stats.txt](https://drive.google.com/file/d/1p-23K2KfFh3_pbZKQj_UCwpWgpD5B-Ew/view?usp=sharing).
Os resultados das filtragem propostas no T1 originaram matrizes de 56 amostras (todos os indivíduos sequenciados) com **1,613,061 SNPs**, sendo que esta matriz de SNPs possui **93,46%** de dados faltantes (missing sites). E matriz de sequências com **35,920,916 nucleotídeos** com **93,46%** de dados faltantes.
###### Resultados de Ugibba_T2:
O **teste T2** foi realizado com a modificação dos seguintes parâmetros: retirar as amostras com 1 indivíduo, sendo assim, das 56 amostras iniciais de *U. gibba* retiramos as amostras 3304, 3368 e 3382; e cada população deve indicar pelo menos 2 *loci* em comum entre os indivíduos para ser considerado válido
Os resultados das filtragens propostas por T2 originaram matrizes de 53 amostras com **9,988 SNPs**, sendo que esta matriz possui **7,37%** de dados faltantes. E a matriz de sequências com **423,152 nucleotídeos**, com **6,02%** de dados faltantes
###### Resultados de Ugibba_T3:
O **teste T3** foi realizado com a modificação dos seguintes parâmetros: retirar as amostras com 2 indivíduos ou menos, sendo assim, das 56 amostras iniciais retiramos as amostras 3304, 3368, 3382, 3307 e 3334; e cada população deve indicar pelo menos 2 *loci* em comum entre os indivíduos para ser considerado válido.
Os resultados das filtragens propostas por T3 originaram matrizes de 49 amostras com **10,617 SNPs**, sendo que esta matriz possui **8,22%** de dados faltantes. E a matriz de sequências com **452,800 nucleotídeos**, com **6,91%** de dados faltantes
###### Resultados de Ugibba_T4:
Os testes seguintes foram feitos com genoma de referência ([Lan et al., 2017](https://www.pnas.org/doi/10.1073/pnas.1702072114)).
### Para gerar o INPUT do Arlequin
:::info
Estes scripts foram retirados e modificados a partir do sítio: [Banta Labs](https://sites.google.com/site/thebantalab/tutorials)
:::
Para calcularmos as Fst nós utilizamos o pragram Arlequin, contudo, a matriz de SNPs tem que ser formatada de acordo para realizar tal análise, sendo assim nós executamos o script em R no Rstudio:
```R!
#set your working directory!
#install.packages("ape")
library(ape)
#to get arlequin format
# read in data
dat <- read.dna("Ugibba_ref_nuclear.snps2.fas", format = "fasta")
dat.matrix <- read.dna("Ugibba_ref_nuclear.snps2.fas", format = "fasta", as.character = TRUE)
if(file.exists("output.arp")){file.remove("output.arp")}
outfile <- cbind(1,rownames(dat.matrix))
colnames(outfile) <- c("group", "name")
write.csv(outfile, "companion data initial.csv", row.names = FALSE)
dat.matrix[dat.matrix == "-"] <- "?"
dat.matrix[dat.matrix == "N"] <- "?"
dat.matrix[dat.matrix == "n"] <- "?"
x2 <- dat.matrix
colnames(x2) <- paste("locus", seq(ncol(x2)), sep = "")
x3 <- rbind(x2, x2)
row.names(x3) <- NULL
x3 <- cbind(x3, 1)
colnames(x3)[ncol(x3)] <- "pop"
for(i in 1:nrow(x2)){
x3[((2*i)-1), 1:(ncol(x3)-1)] <- x2[i,]
x3[((2*i)), 1:(ncol(x3)-1)] <- x2[i, ]
x3[((2*i)-1),"pop"] <- rownames(x2)[i]
x3[((2*i)),"pop"] <- rownames(x2)[i]
}
x4 <- as.matrix(x3)
is.odd <- function(x) x %% 2 != 0
is.even <- function(x) x %% 2 == 0
for(i in 1:(ncol(x4)-1)){
for(j in 1:nrow(x4)){
if(!is.na(x4[j,i])){
if(x4[j,i] == "y"){
if(is.odd(j)){
x4[j,i] <- "c"
}else{
x4[j,i] <- "t"
}
}
if(x4[j,i] == "r"){
if(is.odd(j)){
x4[j,i] <- "a"
}else{
x4[j,i] <- "g"
}
}
if(x4[j,i] == "w"){
if(is.odd(j)){
x4[j,i] <- "a"
}else{
x4[j,i] <- "t"
}
}
if(x4[j,i] == "s"){
if(is.odd(j)){
x4[j,i] <- "c"
}else{
x4[j,i] <- "g"
}
}
if(x4[j,i] == "k"){
if(is.odd(j)){
x4[j,i] <- "t"
}else{
x4[j,i] <- "g"
}
}
if(x4[j,i] == "m"){
if(is.odd(j)){
x4[j,i] <- "c"
}else{
x4[j,i] <- "a"
}
}
}
}
}
rownames(x4) <- x4[,"pop"]
x4 <- x4[order(rownames(x4)),]
x5 <- x4[,-which(colnames(x4) == "pop")]
x6 <- NULL
for(i in 1:nrow(x5)){
x6 <- rbind(x6, noquote((paste(x5[i,], collapse = " "))))
}
names(x6) <- rownames(x5)
if(file.exists("companion data.csv")){cdata <- read.csv("companion data.csv")}else{cdata <- read.csv("companion data initial.csv")}
crlf <- "\r\n"
outfile = "output.arp"
cat('[Profile]', file=outfile, append=FALSE, crlf)
cat('', file=outfile, append=TRUE, crlf)
cat('Title="data"', file=outfile, append=TRUE, crlf)
write.table(paste("NBSamples=",length(table(cdata[,1])), sep = ""), append = TRUE, "output.arp", row.names = FALSE, col.names = FALSE, quote = FALSE, crlf)
cat("", file=outfile, append=TRUE, crlf)
cat('DataType=DNA', file=outfile, append=TRUE, crlf)
cat('GenotypicData=1', file=outfile, append=TRUE, crlf)
cat('GameticPhase=0', file=outfile, append=TRUE, crlf)
cat('LocusSeparator=WHITESPACE', file=outfile, append=TRUE, crlf)
cat("", file=outfile, append=TRUE, crlf)
cat("[Data]", file=outfile, append=TRUE, crlf)
cat("[[Samples]]", file=outfile, append=TRUE, sep = "")
for(i in 1:length(names(table(cdata[,1])))){
cat("", file=outfile, append=TRUE, crlf)
cat("", file=outfile, append=TRUE, crlf)
cat('SampleName=', file=outfile, append=TRUE, sep = "")
write.table(names(table(cdata[,1]))[i], append = TRUE, "output.arp", row.names = FALSE, col.names = FALSE, quote = TRUE, sep = "", crlf)
write.table(paste('SampleSize=', length(which(tolower(cdata[,1]) == tolower(names(table(cdata[,1])))[i])), sep = ""), append = TRUE, "output.arp", row.names = FALSE, col.names = FALSE, quote = FALSE, sep = "", crlf)
cat('SampleData={', file=outfile, append=TRUE, sep = "")
cat("", file=outfile, append=TRUE, crlf)
cdata.sub <- cdata[which(tolower(cdata[,1]) == tolower(names(table(cdata[,1]))[i])),]
for(j in 1:nrow(cdata.sub)){
cat("", file=outfile, append=TRUE, crlf)
x6.sub <- x6[which(names(x6) == cdata.sub[j,2])]
write.table(paste(tolower(cdata.sub[j,2]), 1, toupper(x6.sub[1]), sep = " "), append = TRUE, "output.arp", row.names = FALSE, col.names = FALSE, quote = FALSE, sep = "", crlf)
write.table(toupper(x6.sub[2]), append = TRUE, "output.arp", row.names = FALSE, col.names = FALSE, quote = FALSE, sep = "", crlf)
}
cat("", file=outfile, append=TRUE, crlf)
cat('}', file=outfile, append=TRUE, sep = "")
}
cat('', file=outfile, append=TRUE, crlf)
cat('', file=outfile, append=TRUE, crlf)
cat('[[Structure]]', file=outfile, append=TRUE, crlf)
cat('StructureName="Test"', file=outfile, append=TRUE, crlf)
cat('NbGroups=1', file=outfile, append=TRUE, crlf)
cat('', file=outfile, append=TRUE, crlf)
cat('Group={', file=outfile, append=TRUE, crlf)
cat('', file=outfile, append=TRUE, crlf)
for(i in 1:length(names(table(cdata[,1])))){
write.table(names(table(cdata[,1]))[i], append = TRUE, "output.arp", row.names = FALSE, col.names = FALSE, quote = TRUE, sep = "", crlf)
}
cat('', file=outfile, append=TRUE, crlf)
cat('}', file=outfile, append=TRUE, sep = "")
```
Após executar nós editamos o arquivo **companion data initial.csv** para atribuirmos as populações aos indivíduos e salvamos o arquivo com o nome **"companion data.csv"** e executamos novamente o mesmo script anterior.
Que vai produzir um arquivo **.arp**. Este arquivo é o arquivo de entrada para o Arlequin.
### Cálculo do FST pairwise
Para calcularmos o FST nós carregamos o arquivo de entrada no Arlequin e usamos as opções Settings -> Population comparisons.
Mudamos a opção "Nos of permutations" para aumentarmos o número de permutas da análise para 10000:
E apertamos **"Start"**
### Gráfico para visualizar os FSTs
Para a fácil visualização dos resultados par-a-par nós geramos um HeatMap a partir do script R:
```R!=
#This script comes from: https://nicolawongwaiyee.wordpress.com/2017/01/24/extract-fst-mat-data-from-arlequin-xml-output/
#pkg_list = c("XML","corrplot","magrittr","dplyr")
#check_pkg <- function(pkg){
# if(!require(pkg)){
# install.packages(pkg)
# }
#}
#lapply(pkg_list,check_pkg)
#Load the packages
lapply(pkg_list,require,character.only=TRUE)
arqxml <- xmlParse("output.xml")
arq_fstmat <- xpathSApply(arqxml,"//PairFstMat",xmlValue) %>% as.vector()
arq_fstmat <-strsplit(arq_fstmat,"\n")
arq_fstmat_mt <- do.call(cbind,arq_fstmat[1])
arq_fstmat_mt %<>% subset(arq_fstmat_mt[,1] != "")
arq_fstmat_mt <- gsub(" + "," ",arq_fstmat_mt)
arq_fstmat_mt <- strsplit(arq_fstmat_mt," ")
try(arq_fstmat_mt <- do.call(rbind,arq_fstmat_mt))
arq_fstmat_mt <- arq_fstmat_mt[3:nrow(arq_fstmat_mt),3:ncol(arq_fstmat_mt)]
#Full matrix
arq_fstmat_mt[upper.tri((arq_fstmat_mt))] <- t(arq_fstmat_mt)[upper.tri(arq_fstmat_mt)]
diag(arq_fstmat_mt) <-0
arq_rowname <- xpathSApply(arqxml,"//pairDistPopLabels",xmlValue) %>%
strsplit("\n")
arq_rowname <- do.call(cbind,arq_rowname)
arq_rowname <- gsub("+\\d+:\t","",arq_rowname)
arq_rowname <- arq_rowname[5:nrow(arq_rowname)] %>% as.vector()
arq_fstmat_mt %<>% apply(c(1,2),as.numeric)
rownames(arq_fstmat_mt) <- arq_rowname
colnames(arq_fstmat_mt) <- arq_rowname
#Now the matrix looks...
arq_fstmat_mt[1:5,1:5]
write.csv(arq_fstmat_mt, "arq_fstmat_mt.csv")
arq_fstmat_mt[arq_fstmat_mt<0] <- 0
fstmat_plot <- corrplot(arq_fstmat_mt,method="color",type="lower",na.label = " ",
tl.col="black",tl.srt=8,order="FPC",cl.lim = c(0,1))
#for p-values
arq_fstp <- xpathSApply(arqxml,"//PairFstPvalMat",xmlValue) %>% as.vector()
arq_fstp<-strsplit(arq_fstp,"\n")
arq_fstp <- do.call(cbind,arq_fstp[1])
arq_fstp %<>% subset(arq_fstp[,1] != "")
arq_fstp <- gsub(" + "," ",arq_fstp)
arq_fstp <- strsplit(arq_fstp," ")
try(arq_fstp <- do.call(rbind,arq_fstp))
arq_fstp <- arq_fstp[2:nrow(arq_fstp),3:ncol(arq_fstp)]
arq_fstp[upper.tri(arq_fstp)] <-NA
diag(arq_fstp) <- 0
#Remove the SD of pvalue
arq_fstp <- gsub(" ","",arq_fstp)
arq_fstp <- gsub("\\+.*","",arq_fstp)
write.csv(arq_fstp, "arq_fstp.csv")
#Full matrix
arq_fstp[upper.tri(arq_fstp)] <- t(arq_fstp)[upper.tri(arq_fstp)]
rownames(arq_fstp) <- arq_rowname
colnames(arq_fstp) <- arq_rowname
arq_fstp %<>% apply(c(1,2),as.numeric)
#FstMat with insig p-value mark
corrplot(arq_fstmat_mt,method="color",type="lower",na.label = " ",
tl.col="black",tl.srt=10,order="FPC",cl.lim = c(0,1),
p.mat = arq_fstp,sig.level = 0.05,insig="pch",pch.cex = 1,pch.col = "#6b717a")
```
:::info
As vezes o Arlequin produz esse resultado junto com o correto:
```
<A NAME="27_01_23at17_28_21_gr_pairw_diff"></A>
<data>
======================================================================================
== Comparisons of pairs of population samples
======================================================================================
List of labels for population samples used below:
-------------------------------------------------
</data> <pairDistPopLabels time="27/01/23 at 17:28:21">
Label Population name
----- ---------------
1: Group1
</pairDistPopLabels> <data>
------------------------
Population pairwise FSTs
------------------------
</data> <PairFstMat time="27/01/23 at 17:28:21">
Distance method: Pairwise difference
1
1 0.00000
</PairFstMat> <data>
------------
FST P values
------------
Number of permutations : 10100
</data> <PairFstPvalMat time="27/01/23 at 17:28:21">
1
1 *
</PairFstPvalMat> <data>
------------
Matrix of significant Fst P values
Significance Level=0.05000
------------
Number of permutations : 10100
1
1
</data>
```
Caso acontecer, retirar esse pedaço do arquivo XML.
:::
Sign in with Wallet
Connect another wallet