---
title: "07-CG-motifs"
author: "Kathleen Durkin"
date: "2025-05-16"
always_allow_html: true
output: 
  github_document:
    toc: true
    toc_depth: 3
    number_sections: true
    html_preview: true 
  bookdown::html_document2:
    theme: cosmo
    toc: true
    toc_float: true
    number_sections: true
    code_folding: show
    code_download: true
---

[Assignment details](https://sr320.github.io/course-fish546-2025/assignments/07-CG.html)

Read in one of my fastas (if necessary)
```{r, engine='bash', eval=FALSE}
curl -L https://gannet.fish.washington.edu/gitrepos/ceasmallr/output/00.00-trimming-fastp/CF02-CM02-Zygote_R2_001.fastp-trim.fq.gz -o ../../project/data/trimmed-WGBS-reads/CF02-CM02-Zygote_R2_001.fastp-trim.fq.gz
```

Unzip fasta file and select a subset (the full file is prohibitively large)
```{r, engine='bash', eval=FALSE}
gunzip ../../project/data/trimmed-WGBS-reads/CF02-CM02-Zygote_R2_001.fastp-trim.fq.gz

head -n 4000 ../../project/data/trimmed-WGBS-reads/CF02-CM02-Zygote_R2_001.fastp-trim.fq.gz > ../data/07-CG-motifs/input.fq
```

Load R packages
```{r}
library(seqinr)
library(ShortRead)
```


```{r, cache=TRUE}
# Replace 'input.fasta' with the name of your multi-sequence fasta file
input_file <- "../data/07-CG-motifs/input.fq"

# Need to use silghtly different code to accomodate an .fq file instead of a .fasta
# sequences <- read.fasta(input_file)
fq <- readFastq(input_file)
sequences <- sread(fq)

```


```{r}
# Set the seed for reproducibility (optional)
set.seed(42)

number_of_sequences_to_select <- 100

if (length(sequences) < number_of_sequences_to_select) {
  warning("There are fewer than 10 sequences in the fasta file. All sequences will be selected.")
  number_of_sequences_to_select <- length(sequences)
}

selected_indices <- sample(length(sequences), number_of_sequences_to_select)
selected_sequences <- sequences[selected_indices]

```


```{r}
# Replace 'output.fasta' with your desired output file name
output_file <- "../output/07-CG-motifs/100_seqs.fasta"
write.fasta(selected_sequences, names(selected_sequences), output_file, open = "w")
```

Issues with IGV visualization because this generated fasta doesn't have matching scaffold name. Fix this here
```{r, engine='bash'}
sed -i '1s/.*/>EMBOSS_001/' ../output/07-CG-motifs/100_seqs.fasta
```


```{bash}
#needed downstream for IGV
/home/shared/samtools-1.12/samtools faidx \
../output/07-CG-motifs/100_seqs.fasta
```


```{bash}
fuzznuc -sequence ../output/07-CG-motifs/100_seqs.fasta -pattern CG -rformat gff -outfile ../output/07-CG-motifs/CGoutput.gff
```


IGV visualization:

![](../output/07-CG-motifs/IGV.png)


