Multiple sequence alignment Fasta file manipulation

# (1) Read in sample information
source_mods = readxl::read_xlsx("G:/My Drive/Tetramesa Project/Source Modifiers/source_modifiers_updated_10_06_2021.xlsx", 
                                sheet = "cleaned_tetramesa_clarke"
)

# (2) Read in binary values
host.plants = readxl::read_xlsx("G:/My Drive/Tetramesa Project/Source Modifiers/source_modifiers_updated_10_06_2021.xlsx", 
                                sheet = "host_plants"
)

# Change any NA cells to contain zero values
host.plants[is.na(host.plants)] = 0

# Change the values in the dataframe to become characters (in case they are read as strings)
host.plants[] = lapply(host.plants, as.character)

# If the file is read into R using the readxl:: readxlsx() function, then this line changes the first column to rownames
host.plants = tibble::column_to_rownames(host.plants, var = "sample_id")

# Clean up the column names
host.plants = janitor::clean_names(host.plants)

# Subset the sample information datasheet accordingly. In my data, the columns after spyramidalis contain other information, which isn't needed here. 
# I only want to look at the Arundo donax (adonax) to Sporobolus pyramidalis (spyramidalis) columns
host.plants.filter = dplyr::select(host.plants, adonax:spyramidalis)

# (3) Read in the MSA fasta file
seqs_28S = ape::read.FASTA("C:/Users/s1000334/Documents/PhD_Tetramesa/Raw Sequences/Tetramesa/ALIGNED_READY/Tetramesa/Single genes/More_eurytomidae_included/28S/28S_only_aligned.fasta")

##########################
# FILTERING
##########################

# I want to subset sequences for only those that belong to the Golden shouldered (GS) morphotype from my sample info sheet (source_mods):
GS = dplyr::filter(source_mods, morphotype == "GS")

# Create an empty fasta file to which the subsetted sequences will be written
file.create("28S_GS.fasta")

# Loop through each sequence_ID name in the subsetted GS dataframe, and extract the corresponding sequence name from the uploaded Fasta alignment file.
# Then write that sequence to a new file. In this way, you will end up with a subsetted Fasta alignment with only the golden shouldered specimens

for(i in 1:nrow(GS)){
  
  for(j in 1:length(seqs_28S)){
    if(GS$sequence_ID[i] == labels(seqs_28S[j]))
      ape::write.FASTA( seqs_28S[j], "28S_GS.fasta", append = T )
  }
  
}

###############################################################
# CREATE A TRAITS BLOCK FOR POPART (HAPLOTYPE NETWORK PROGRAM)
###############################################################

# Using the newly subsetted Fasta alignment, we now want to subset the binary values matrix (host.plants) to match it

# Read in the subsetted Fasta file created above
seqs_28S_GS = ape::read.FASTA("C:/Users/s1000334/Documents/PhD_Tetramesa/Raw Sequences/Tetramesa/ALIGNED_READY/Tetramesa/Single genes/More_eurytomidae_included/28S/PopArt/28S_GS.fasta")

# Extract the rows in the binary data frame (host.plants) for the only the sample names that are present in the subsetted Fasta file (seqs_28S_GS)
host.plants.popart.block = host.plants.filter %>% dplyr::filter(row.names(host.plants.filter) %in% labels(seqs_28S_GS))

# Get the binary data frame into the same order as the sequence names in the subsetted Fasta file
host.plants.popart.block.ordered = host.plants.popart.block[match(labels(seqs_28S_GS), rownames(host.plants.popart.block)),]

# Find which columns have only zeros, so that those can be removed (no point in keeping a trait column if there are no sequence representatives in it!)
zero_cols = c()
for(i in 1:ncol(host.plants.popart.block.ordered)){
  if( all(host.plants.popart.block.ordered[i] == 0) )
    zero_cols[i] = i
}

# Remove the NA values from the vector
zero_cols = zero_cols[!is.na(zero_cols)]

# Remove the columns that had only zeroes
host.plants.popart.block.ordered = dplyr::select(host.plants.popart.block.ordered, -zero_cols )

# Write the resulting binary matrix to a text file, which can be appended to the end of the sequence alignment Nexus file as a TRAITS block
# Just remove the comma after each sample name to get it into the correct format for PopArt
write.table(host.plants.popart.block.ordered, file = "popart_traits_block.txt", quote = FALSE, col.names = TRUE, sep = ",")