Population label key: The code labels used in the analysis below correspond to the following manuscript labels: AUT = AUTO, MAN = NON-AUTO, NEW = NON-AUTO-FIELD.
Note: The analyses in this file were not used directly in the final version of the manuscript. However, this notebook generates intermediate data files (snps_genes_chr.rds, SNPs_on_DE_genes.txt, SNPs_79_DE.txt) that are consumed by the Linkage Network Analysis (File S7) and Fst (File S8) notebooks.
## Warning in system("timedatectl", intern = TRUE): running command 'timedatectl'
## had status 1 and error message 'Function not implemented'Since we created a chromosomal scale and the gff file has the scaffolds positions, we can get the genomic coordinates of the genes and then lift them to the chromosomal scale. Another way would be get the SNPs IDs and check where they are located. Each way has its up and down sides.
The gff file with the annotation for the AalbF3 genome assembly is available in the project, under the files directory.
## chr3.9 Gnomon gene 10376 142706 . - . ID=gene-LOC109400916;Dbxref=GeneID:109400916;Name=LOC109400916;gbkey=Gene;gene=LOC109400916;gene_biotype=protein_coding;partial=true;start_range=.,10376
## chr3.9 Gnomon mRNA 10376 142706 . - . ID=rna-XM_029851713.1;Parent=gene-LOC109400916;Dbxref=GeneID:109400916,Genbank:XM_029851713.1;Name=XM_029851713.1;gbkey=mRNA;gene=LOC109400916;model_evidence=Supporting evidence includes similarity to: 1 Protein%2C and 97%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 57 samples with support for all annotated introns;partial=true;product=UPF0565 protein C2orf69 homolog;start_range=.,10376;transcript_id=XM_029851713.1
## chr3.9 Gnomon exon 142438 142706 . - . ID=exon-XM_029851713.1-1;Parent=rna-XM_029851713.1;Dbxref=GeneID:109400916,Genbank:XM_029851713.1;gbkey=mRNA;gene=LOC109400916;partial=true;product=UPF0565 protein C2orf69 homolog;transcript_id=XM_029851713.1
## chr3.9 Gnomon exon 60465 60496 . - . ID=exon-XM_029851713.1-2;Parent=rna-XM_029851713.1;Dbxref=GeneID:109400916,Genbank:XM_029851713.1;gbkey=mRNA;gene=LOC109400916;partial=true;product=UPF0565 protein C2orf69 homolog;transcript_id=XM_029851713.1
## chr3.9 Gnomon exon 11273 11676 . - . ID=exon-XM_029851713.1-3;Parent=rna-XM_029851713.1;Dbxref=GeneID:109400916,Genbank:XM_029851713.1;gbkey=mRNA;gene=LOC109400916;partial=true;product=UPF0565 protein C2orf69 homolog;transcript_id=XM_029851713.1
## chr3.9 Gnomon exon 10638 10870 . - . ID=exon-XM_029851713.1-4;Parent=rna-XM_029851713.1;Dbxref=GeneID:109400916,Genbank:XM_029851713.1;gbkey=mRNA;gene=LOC109400916;partial=true;product=UPF0565 protein C2orf69 homolog;transcript_id=XM_029851713.1
## chr3.9 Gnomon exon 10376 10570 . - . ID=exon-XM_029851713.1-5;Parent=rna-XM_029851713.1;Dbxref=GeneID:109400916,Genbank:XM_029851713.1;gbkey=mRNA;gene=LOC109400916;partial=true;product=UPF0565 protein C2orf69 homolog;start_range=.,10376;transcript_id=XM_029851713.1
## chr3.9 Gnomon CDS 11273 11629 . - 0 ID=cds-XP_029707573.1;Parent=rna-XM_029851713.1;Dbxref=GeneID:109400916,Genbank:XP_029707573.1;Name=XP_029707573.1;gbkey=CDS;gene=LOC109400916;partial=true;product=UPF0565 protein C2orf69 homolog;protein_id=XP_029707573.1
## chr3.9 Gnomon CDS 10638 10870 . - 0 ID=cds-XP_029707573.1;Parent=rna-XM_029851713.1;Dbxref=GeneID:109400916,Genbank:XP_029707573.1;Name=XP_029707573.1;gbkey=CDS;gene=LOC109400916;partial=true;product=UPF0565 protein C2orf69 homolog;protein_id=XP_029707573.1
## chr3.9 Gnomon CDS 10376 10570 . - 1 ID=cds-XP_029707573.1;Parent=rna-XM_029851713.1;Dbxref=GeneID:109400916,Genbank:XP_029707573.1;Name=XP_029707573.1;gbkey=CDS;gene=LOC109400916;partial=true;product=UPF0565 protein C2orf69 homolog;protein_id=XP_029707573.1;start_range=.,10376To start we can get the scaffold, start and end of each gene. I will do this because I might need to represent the genes in a plot using ggplot.
# In the code below we get the gene id, the scaffold id, the gene start, the gene end, and the size of the genomic regions the gene covers (end - start)
awk 'BEGIN {OFS="\t"} $3=="gene" {split($9, a, ";"); split(a[1], b, "="); size = $5 - $4; print b[2], $1, $4, $5, size}' data/files/genes.gff | head## gene-LOC109400916 chr3.9 10376 142706 132330
## gene-Trnav-aac chr3.9 33853 33925 72
## gene-Trnat-agu chr3.9 34665 34738 73
## gene-Trnah-gug chr3.9 34864 34935 71
## gene-Trnah-gug-2 chr3.9 35682 35753 71
## gene-Trnah-gug-3 chr3.9 49338 49409 71
## gene-LOC115261967 chr3.9 59372 68451 9079
## gene-Trnav-aac-2 chr3.9 69409 69481 72
## gene-Trnat-agu-2 chr3.9 145258 145328 70
## gene-LOC115254308 chr3.9 206669 220356 13687Create a file with the genomic coordinates for each gene
# In the code below we get the gene id, the scaffold id, the gene start, the gene end, and the size of the genomic regions the gene covers (end - start)
awk 'BEGIN {OFS="\t"} $3=="gene" {split($9, a, ";"); split(a[1], b, "="); size = $5 - $4; print b[2], $1, $4, $5, size}' data/files/genes.gff > data/files/genes.txt;
wc -l data/files/genes.txt## 20621 data/files/genes.txtWe got 20,621 genes in the gff file.
We can import the data into R now
gene_data <- read_delim(here("data", "files","genes.txt"), delim = "\t", col_names = c("Gene_ID", "Scaffold", "Start", "End", "Size"), show_col_types = FALSE)
head(gene_data)## # A tibble: 6 x 5
## Gene_ID Scaffold Start End Size
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 gene-LOC109400916 chr3.9 10376 142706 132330
## 2 gene-Trnav-aac chr3.9 33853 33925 72
## 3 gene-Trnat-agu chr3.9 34665 34738 73
## 4 gene-Trnah-gug chr3.9 34864 34935 71
## 5 gene-Trnah-gug-2 chr3.9 35682 35753 71
## 6 gene-Trnah-gug-3 chr3.9 49338 49409 71First we need to remove the string “chr” from the collumn Scaffold.
## # A tibble: 6 x 5
## Gene_ID Scaffold Start End Size
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 gene-LOC109400916 3.9 10376 142706 132330
## 2 gene-Trnav-aac 3.9 33853 33925 72
## 3 gene-Trnat-agu 3.9 34665 34738 73
## 4 gene-Trnah-gug 3.9 34864 34935 71
## 5 gene-Trnah-gug-2 3.9 35682 35753 71
## 6 gene-Trnah-gug-3 3.9 49338 49409 71Now we can import the bim file in which we have the Scaffold
# Import the function
source(here("notebooks", "helpers", "import_bim.R"))
# Import the data
scaffolds <- import_bim(here("output", "quality_control", "file1.bim"))
head(scaffolds)## # A tibble: 6 x 6
## Scaffold SNP Cm Position Allele1 Allele2
## <chr> <chr> <int> <dbl> <chr> <chr>
## 1 1.1 AX-581444870 0 97856 C T
## 2 1.1 AX-583033226 0 161729 A G
## 3 1.1 AX-583035067 0 229640 T A
## 4 1.1 AX-583035083 0 305518 A G
## 5 1.1 AX-583035102 0 308124 A G
## 6 1.1 AX-583033342 0 315059 C GNow import the file with chromosomal scale
chromosomes <- import_bim(here("output", "quality_control", "file7.bim"))
# Rename columns
chromosomes <- chromosomes |>
dplyr::rename(Chromosome = Scaffold,
Position_chr = Position)
head(chromosomes)## # A tibble: 6 x 6
## Chromosome SNP Cm Position_chr Allele1 Allele2
## <chr> <chr> <int> <dbl> <chr> <chr>
## 1 1 AX-581444870 0 97856 C T
## 2 1 AX-583035083 0 305518 A G
## 3 1 AX-583035102 0 308124 A G
## 4 1 AX-583033342 0 315059 C G
## 5 1 AX-583035163 0 315386 A G
## 6 1 AX-583033356 0 315674 C TNow merge the object
# Merge and keep only the rows with data
snps_albo <- inner_join(scaffolds, chromosomes, by = "SNP", suffix = c("", "")) |>
dplyr::select(
SNP, Chromosome, Position_chr, Scaffold, Position
)
head(snps_albo)## # A tibble: 6 x 5
## SNP Chromosome Position_chr Scaffold Position
## <chr> <chr> <dbl> <chr> <dbl>
## 1 AX-581444870 1 97856 1.1 97856
## 2 AX-583035083 1 305518 1.1 305518
## 3 AX-583035102 1 308124 1.1 308124
## 4 AX-583033342 1 315059 1.1 315059
## 5 AX-583035163 1 315386 1.1 315386
## 6 AX-583033356 1 315674 1.1 315674Now we can add our gene data with the scaffold data
# Perform a non-equi join
merged_scaff <- scaffolds |>
left_join(gene_data, by = "Scaffold") |>
filter(Position >= Start & Position <= End)
head(merged_scaff)## # A tibble: 6 x 10
## Scaffold SNP Cm Position Allele1 Allele2 Gene_ID Start End Size
## <chr> <chr> <int> <dbl> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 1.1 AX-58144~ 0 97856 C T gene-L~ 64508 102530 38022
## 2 1.1 AX-58303~ 0 315059 C G gene-L~ 310594 316608 6014
## 3 1.1 AX-58303~ 0 315386 A G gene-L~ 310594 316608 6014
## 4 1.1 AX-58303~ 0 315674 C T gene-L~ 310594 316608 6014
## 5 1.1 AX-58303~ 0 346460 A T gene-L~ 346347 346851 504
## 6 1.10 AX-58303~ 0 91677 T C gene-L~ 66918 192829 125911Now we merge with chromosomal data
# Merge
snps_genes_chr <- inner_join(merged_scaff, chromosomes, by = "SNP", suffix = c("", "")) |>
dplyr::select(
SNP, Chromosome, Position_chr, Scaffold, Position, Gene_ID, Start, End, End
) |>
mutate(Gene_ID = str_remove(Gene_ID, "gene-"))
head(snps_genes_chr)## # A tibble: 6 x 8
## SNP Chromosome Position_chr Scaffold Position Gene_ID Start End
## <chr> <chr> <dbl> <chr> <dbl> <chr> <dbl> <dbl>
## 1 AX-581444870 1 97856 1.1 97856 LOC10939~ 64508 102530
## 2 AX-583033342 1 315059 1.1 315059 LOC10943~ 310594 316608
## 3 AX-583035163 1 315386 1.1 315386 LOC10943~ 310594 316608
## 4 AX-583033356 1 315674 1.1 315674 LOC10943~ 310594 316608
## 5 AX-583035257 1 442875 1.10 91677 LOC10940~ 66918 192829
## 6 AX-583035268 1 462944 1.10 111746 LOC10940~ 66918 192829Save it
Now we can import gene expression data
gene_expression <- read_delim(here("data", "files","MANvsAUTO_sig_mRNAs.csv"), delim = ",", col_names = TRUE, show_col_types = FALSE) |>
dplyr::select(
gene,log2FoldChange
) |>
dplyr::rename(
Gene_ID = gene
)
head(gene_expression)## # A tibble: 6 x 2
## Gene_ID log2FoldChange
## <chr> <dbl>
## 1 LOC115262812 -10.5
## 2 LOC109401291 -9.19
## 3 LOC109397830 8.73
## 4 LOC115264022 8.56
## 5 LOC115260314 8.51
## 6 LOC115258723 8.39Now we can merge it with our snps_genes_chr object
snps_expression <- inner_join(snps_genes_chr, gene_expression, by = "Gene_ID", suffix = c("", ""))
head(snps_expression)## # A tibble: 6 x 9
## SNP Chromosome Position_chr Scaffold Position Gene_ID Start End
## <chr> <chr> <dbl> <chr> <dbl> <chr> <dbl> <dbl>
## 1 AX-583050970 1 3777633 1.10 3426435 LOC10942~ 3.42e6 3.43e6
## 2 AX-583052853 1 3777864 1.10 3426666 LOC10942~ 3.42e6 3.43e6
## 3 AX-583051011 1 3779183 1.10 3427985 LOC10942~ 3.42e6 3.43e6
## 4 AX-583052887 1 3793314 1.10 3442116 LOC10942~ 3.44e6 3.47e6
## 5 AX-583051025 1 3793653 1.10 3442455 LOC10942~ 3.44e6 3.47e6
## 6 AX-583052920 1 3794185 1.10 3442987 LOC10942~ 3.44e6 3.47e6
## # i 1 more variable: log2FoldChange <dbl>Check how many SNPs we have in the genes with differential expression
## [1] 1192Save the list of SNPs on DE genes
# Select only the desired columns
data_to_save <- snps_expression[, c("SNP", "Chromosome", "Position_chr")]
# Write the selected data to a file
write.table(
data_to_save,
file = here("output", "snpeff", "SNPs_on_DE_genes.txt"),
row.names = FALSE,
quote = FALSE,
col.names = TRUE,
sep = "\t" # or use "," for CSV format
)Now we can check if any of the SNPs we found in our analysis are within these genes. If we do not find any, them we can develop a different approach, which is find the nearest gene for the SNPs or the nearest SNP for the genes of interest.
First lets check the 17 SNPs that we find in common
snps01<- read_delim(here("output", "pcadapt", "4_way_venn_common_SNPs_pcadapt_outflank.txt"), delim = "\t", col_names = c("SNP"), show_col_types = FALSE)
head(snps01)## # A tibble: 6 x 1
## SNP
## <chr>
## 1 AX-583515734
## 2 AX-585196879
## 3 AX-579548089
## 4 AX-579560686
## 5 AX-579584369
## 6 AX-579583040Check how many SNPs
## [1] 16Check if there is any differentially expressed gene with one of these SNPs
# Filter snps_expression to keep only the rows that have a matching SNP in snps01
filtered_snps_expression <- snps_expression |>
semi_join(snps01, by = "SNP")
head(filtered_snps_expression)## # A tibble: 0 x 9
## # i 9 variables: SNP <chr>, Chromosome <chr>, Position_chr <dbl>,
## # Scaffold <chr>, Position <dbl>, Gene_ID <chr>, Start <dbl>, End <dbl>,
## # log2FoldChange <dbl>None of the 17 SNPs are present in the genes. Now we can expand our search.
# Lets check all the SNPs from the comparisons of the populations vs the autogenous line
# Read data from txt files
MAN_AUT <-
read.table(
here("output", "outflank", "man_aut_SNPs_outFlank.txt"),
stringsAsFactors = FALSE
) |>
drop_na()
NEW_AUT <-
read.table(
here("output", "outflank", "new_aut_SNPs_outFlank.txt"),
stringsAsFactors = FALSE
)
NEW_MAN_AUT <-
read.table(
here("output", "outflank", "SNPs_outFlank.txt"),
stringsAsFactors = FALSE
)Merge them and get a SNP count
# Combine all three data frames into one
combined_df <- bind_rows(MAN_AUT, NEW_AUT, NEW_MAN_AUT)
# Count unique values and add a 'Count' column
result <- combined_df %>%
group_by(V1) %>%
summarise(Count = n()) %>%
arrange(V1) |>
dplyr::rename(SNP = V1)
# Show the result
head(result)## # A tibble: 6 x 2
## SNP Count
## <chr> <int>
## 1 AX-579436928 1
## 2 AX-579437189 1
## 3 AX-579437701 1
## 4 AX-579438128 2
## 5 AX-579438621 2
## 6 AX-579438851 1How many SNPs?
## [1] 8018How many SNPs out of the 158 are on genes with differential expression
filtered_snps_expression2 <- snps_expression |>
semi_join(result, by = "SNP")
length(filtered_snps_expression2$SNP)## [1] 76## [1] "AX-583054970" "AX-583055828" "AX-583093532" "AX-583142560" "AX-583237406"
## [6] "AX-583279927" "AX-583320377" "AX-583431590" "AX-583462099" "AX-583683199"
## [11] "AX-583687374" "AX-583696963" "AX-583709904" "AX-583778551" "AX-583972796"
## [16] "AX-584188695" "AX-584222538" "AX-584275318" "AX-584506677" "AX-585047348"
## [21] "AX-584538667" "AX-584629321" "AX-584653374" "AX-584653471" "AX-585162248"
## [26] "AX-585162361" "AX-585167663" "AX-584659072" "AX-584659195" "AX-585191554"
## [31] "AX-582458673" "AX-582481773" "AX-582481882" "AX-582486287" "AX-582580299"
## [36] "AX-584792883" "AX-584813246" "AX-585454571" "AX-584949353" "AX-579485773"
## [41] "AX-579518040" "AX-579519594" "AX-579528676" "AX-579530323" "AX-579551827"
## [46] "AX-579587308" "AX-579602820" "AX-580037979" "AX-580363853" "AX-580385664"
## [51] "AX-580546704" "AX-580706309" "AX-580713461" "AX-580734190" "AX-580790145"
## [56] "AX-580878016" "AX-581014751" "AX-581302901" "AX-581487305" "AX-581504020"
## [61] "AX-581505925" "AX-581506313" "AX-581504582" "AX-581506954" "AX-581798715"
## [66] "AX-581839987" "AX-582039669" "AX-582287268" "AX-582293283" "AX-582375521"
## [71] "AX-582683159" "AX-582740296" "AX-582743046" "AX-582772850" "AX-582812598"
## [76] "AX-582813926"Save it
write.table(
filtered_snps_expression2$SNP,
file = here("output", "snpeff","SNPs_79_DE.txt"),
row.names = FALSE,
quote = FALSE,
col.names = FALSE,
sep = "\n"
)We can also check if any of the 157 outliers
Select the 17 SNPs outliers identified in all 3 comparisons where we had the AUT line. See the 4 way Venn diagram
snps_17 <-
read.table(
here("output", "pcadapt", "4_way_venn_common_SNPs_pcadapt_outflank.txt"),
stringsAsFactors = FALSE
)
# Get the 17 SNPs
snps_aut17 <- snps_expression |>
filter(SNP %in% snps_17$V1)
head(snps_aut17)## # A tibble: 0 x 9
## # i 9 variables: SNP <chr>, Chromosome <chr>, Position_chr <dbl>,
## # Scaffold <chr>, Position <dbl>, Gene_ID <chr>, Start <dbl>, End <dbl>,
## # log2FoldChange <dbl>snps_157 <-
read.table(
here("output", "snpeff", "SNPs_158.txt"),
stringsAsFactors = FALSE
)
# Get the 17 SNPs
snps_aut157 <- snps_expression |>
filter(SNP %in% snps_157$V1)
head(snps_aut157)## # A tibble: 4 x 9
## SNP Chromosome Position_chr Scaffold Position Gene_ID Start End
## <chr> <chr> <dbl> <chr> <dbl> <chr> <dbl> <dbl>
## 1 AX-583054970 1 5342242 1.10 4991044 LOC10942~ 4.91e6 5.02e6
## 2 AX-583972796 1 301386519 1.8 315751 LOC10940~ 3.14e5 3.16e5
## 3 AX-581302901 3 191256607 3.170 12036013 LOC10940~ 1.20e7 1.20e7
## 4 AX-581504582 3 226420377 3.19 1649039 LOC10940~ 1.33e6 1.73e6
## # i 1 more variable: log2FoldChange <dbl># Perform an inner join using the SNP column from snps_genes_chr and the V1 column from snps_157
snps_157_genes <- inner_join(snps_genes_chr, snps_157, by = c("SNP" = "V1"))
# View the first few rows of the merged dataset
head(snps_157_genes)## # A tibble: 6 x 8
## SNP Chromosome Position_chr Scaffold Position Gene_ID Start End
## <chr> <chr> <dbl> <chr> <dbl> <chr> <dbl> <dbl>
## 1 AX-583054970 1 5342242 1.10 4991044 LOC10942~ 4.91e6 5.02e6
## 2 AX-583095890 1 19188676 1.105 311276 LOC10941~ 3.04e5 4.55e5
## 3 AX-583355868 1 100679287 1.141 118038 LOC11525~ 3.25e4 1.52e5
## 4 AX-583423493 1 120579795 1.147 4673825 LOC10940~ 4.44e6 4.81e6
## 5 AX-583426050 1 120874430 1.147 4968460 LOC10941~ 4.96e6 4.97e6
## 6 AX-583467551 1 130683717 1.149 2757434 LOC10942~ 2.60e6 2.76e6Save it
write.xlsx(
snps_157_genes,
here(
"output", "snpeff","genes_with_114_SNPs_out_of_157_outliers.xlsx"
)
)Lets look what are the nearest gene of all the 158 SNPs, then we can filter out whatever way we want. We need to do SNP functional annotation for all the SNPs.
To do the functional annotation we need to include all SNPs, so we can use file1.
Let’s create a vcf file for the annotation
plink2 \
--allow-extra-chr \
--bfile output/quality_control/file1 \
--recode vcf \
--out output/snpeff/albo \
--silent;
grep "samples\|variants" output/snpeff/albo.log## 61 samples (25 females, 26 males, 10 ambiguous; 61 founders) loaded from
## 129381 variants loaded from output/quality_control/file1.bim.Now check the R Markdown with functional annotation
Create bed file
snps_functional <- snps_genes_chr %>%
filter(SNP %in% result$SNP) |>
dplyr::select(
SNP, Scaffold, Position
)
head(snps_functional)## # A tibble: 6 x 3
## SNP Scaffold Position
## <chr> <chr> <dbl>
## 1 AX-583035257 1.10 91677
## 2 AX-583035671 1.10 958271
## 3 AX-583037587 1.10 976678
## 4 AX-583039330 1.10 1143681
## 5 AX-583040214 1.10 1462572
## 6 AX-583043295 1.10 1649719# Assuming snps_functional
write.table(snps_functional %>%
select(SNP, Scaffold, Position) %>%
mutate(End = Position + 1),
file = here("output", "snpeff","snps_functional.bed"),
quote = FALSE,
row.names = FALSE,
col.names = FALSE,
sep = "\t")Check it out
## AX-583035257 1.10 91677 91678
## AX-583035671 1.10 958271 958272
## AX-583037587 1.10 976678 976679
## AX-583039330 1.10 1143681 1143682
## AX-583040214 1.10 1462572 1462573
## AX-583043295 1.10 1649719 1649720
## AX-583044718 1.10 1890345 1890346
## AX-583046569 1.10 2585730 2585731
## AX-583051721 1.10 3198232 3198233
## AX-583054970 1.10 4991044 4991045Note: The cluster-gene intersection analysis (clusters 14 and 6 with differentially expressed genes) has been moved to File S4e (LDna combined analysis), which creates the cluster objects and performs the intersection in a single notebook.