3mma-mack
commented
3 years ago First step was to combine the three fasta files for the 5UTRs, cds sequence and 3UTRs into one file which for each gene had the sequences in the order 5UTR-CDS-3UTR. The original plan is shown in the pseudocode below:
# Take the 5UTR, CDS and 3UTR sequences
# Stick them together
# Produce a new fasta file that has sequences in the format 5UTR-CDS-3UTR
# inputs
# S_pombe 5UTR fasta, CDS fasta, 3UTR fasta.
# these should have sequences for individual genes/exons, labelled with IDs
# Downloaded from PomBase
# Process
# check that sequences are labelled with IDs
# for loop to add 5UTR to CDS
# for sequence i in CDS file
# store the ID of sequence i
# find the sequence with the matching ID in 5UTRs
# some may not have 5UTRs, if there isn't then i ++
# merge the sequences into the format 5UTR-CDS
# store product in position i of new object 5UTR-CDS
# return 5UTR-CDS
# i ++
# repeat above for adding 3UTR to 5UTR-CDS
# for sequence i in 5UTR-CDS file
# store the ID of sequence i
# find the sequence with the matching ID in 3UTRs
# some may not have 5UTRs, if there isn't then i ++
# merge the sequences into the format 5UTR-CDS-3UTR
# store product in position i of new object 5UTR-CDS-3UTR
# return 5UTR-CDS-3UTR
# i ++This was adapted and the fasta file was created and checked with the following code:
library(Biostrings)
library(rtracklayer)
library(stringr)
library(GenomicRanges)
library(parallel)
library(rhdf5)
three_UTR <- readDNAStringSet('3UTR.fa.gz')
CDS <- readDNAStringSet('cds.fa.gz')
five_UTR <- readDNAStringSet('5UTR.fa.gz')
cc <- 1
fiveutr_CDS_threeutr <- list()
# loop to add, if present, 5UTRs upstream of CDS sequences and 3UTRs downstream.
# add names column
for (i in unique(names(CDS))){
fiveutr_CDS_threeutr[[cc]] <- DNAStringSet(paste0(unlist(five_UTR[names(five_UTR)==i]),
unlist(CDS[names(CDS)==i]),
unlist(three_UTR[names(three_UTR)==i])
))
cc <- cc + 1
}
fiveutr_CDS_threeutr <- unlist(DNAStringSetList(fiveutr_CDS_threeutr))
names(fiveutr_CDS_threeutr) <- unique(names(CDS))
store <- fiveutr_CDS_threeutr
# check that the length of a sequence fiveUTR_CDS_seqlist is equal
# to the combined lengths of seperate CDS and five_utr length for each gene
for (j in unique(names(fiveutr_CDS_threeutr))){
if (width(fiveutr_CDS_threeutr[names(fiveutr_CDS_threeutr)==j])!=
(sum(width(five_UTR[names(five_UTR)==j])) +
sum(width(CDS[names(CDS)==j])) +
sum(width(three_UTR[names(three_UTR)==j])))){
print('Length of full sequence is not the sum of the length of the 3UTR, CDS and 5UTR')
break
}
}
# names used in gff and fasta have different format, adding '.1' to fasta to make consistent
names(fiveutr_CDS_threeutr) <- paste(names(fiveutr_CDS_threeutr), '.1', sep = '')
# export fasta files with 3UTR, CDS and 5UTR sequences.
writeXStringSet(fiveutr_CDS_threeutr,filepath = file.path('.','S_pombe_full_UTR.fasta'),format = "fasta")
For my project I am wanting to look at how specific genes are regulated, so will likely need gff and fasta files that include full UTRs rather than just a buffer. Pombase has separate fasta files for CDS, 5UTR and 3UTR regions, so I am going to take these and combine them to make new fasta and gff files. Doing this will allow me to look at uOFRs involved in regulation, particularly in genes such as FIL1. However, as the lengths of sequences surrounding the CDS will be inconsistent there may be some issues with the buffer parameter in the config file