Count CG motifs in Roslin genome

[yaaminiv]

I want to count the number of CG motifs in the Roslin genome. I used the following command to do so:

~4 million CpGs seems low (C. virginica has ~14 million). When I look at methylKit::unite output, I get ~4.5 million CpGs with data. Clearly I'm counting the CGs incorrectly. What am I missing?

edited to include unite information`

[kubu4]

This will miss instances when a CG occurs at a newline.

e.g.:

AATCCGATAGC
GATTTACTATCG

[yaaminiv]

Ah! Looked into my old code and am using fgrep -o -i CG {genome-file} | wc -l instead!

[kubu4]

I don't think that will count across newlines, though...

[yaaminiv]

It might? I got ~13 million CG motifs. I will try the tab-delimited approach and compare!

[yaaminiv]

For reference, here's the issue where we had a discussion about fgrep: #685 (comment)

[kubu4]

Ah! Funny we've discussed this before!

[kubu4]

I think the approach is to convert the FastA into a tab-delimited file. I use seqkit to do this, usually. Then, count with awk (possibly need to use gawk, to accommodate large number of fields). E.g.

  # Use seqkit to convert FastA to tab-delimited and print sequence length 
  ${seqkit} fx2tab \
  --length \
  ${fasta} \
  > ${fasta}.tab

  # Print only sequences to new file
  gawk '{ print $2 }' ${fn}_tab > ${fasta}.tab2
  
  # Delimit sequences on CGs and print the number of fields minus 1 to get the number of CGs present.
  gawk -F\[Cc][Gg] '{print NF-1}' ${fasta}.tab2 > CG

[yaaminiv]

Here's the code I'm using (I modified the last chunk to count the number of CGs instead of printing how many there were on each line:

I get 0 CGs. I think this is because the seqkit fx2tab is not working:

How should I modify my code?

[sr320]

You could convert using this instead

#Converting FASTA to tabular format and placing output file in analyses directory
!perl -e '$count=0; $len=0; while(<>) {s/\r?\n//; s/\t/ /g; if (s/^>//) { if ($. != 1) {print "\n"} s/ |$/\t/; $count++; $_ .= "\t";} else {s/ //g; $len += length($_)} print $_;} print "\n"; warn "\nConverted $count FASTA records in $. lines to tabular format\nTotal sequence length: $len\n\n";' \
33496_Ahyacinthus_CoralContigs.fasta > ../../analyses/Ahya/fasta2tab

[yaaminiv]

Still running into an error with conversion to tab...when I look at head it prints the entire file

[yaaminiv]

I might abandon this since I just wanted a rough count of the CGs in the genome. I will make a CG motif track for IGV as well, which will then get me an exact count!

[yaaminiv]

@sr320 How do you usually generate the CG motif files for genomes? Perhaps that's what I should be doing, so I have that as an IGV track

[sr320]

no EMBOSS code - used Galaxy

[yaaminiv]

Does this look right?

[sr320]

Yes (I have it running to now on Galaxy) hopefully we get same product

[sr320]

Here is a link to the GFF https://usegalaxy.org/datasets/bbd44e69cb8906b5bba9b46c41a81f2b/display?to_ext=gff

[kubu4]

Here's fuzznuc code I've used in the past:

https://robertslab.github.io/sams-notebook/2019/08/21/Data-Wrangling-Create-a-CpG-GFF-from-Pgenerosa_v074-using-EMBOSS-fuzznuc-on-Swoose.html

~/programs/EMBOSS-6.6.0/emboss/fuzznuc \
-sequence ~/data/P_generosa/genomes/Pgenerosa_v074.fa \
-pattern CG \
-outfile 20190821_fuzznuc_pgenv074.gff \
-rformat gff

[sr320]

Also note I will separately be dealing with the general ability to count CG at #1207

[kubu4]

normally EMBOSS Fuzznuc....

This is definitely the way to go!

[yaaminiv]

Here's my CG motif track (will add to the handbook): http://owl.fish.washington.edu/halfshell/genomic-databank/cgigas_uk_roslin_v1_fuzznuc_CGmotif.gff Similar count to the 13080574 from fgrep

Looks the same as Steven's: