Manipulating data in R using the Tidyverse

I've just come across the Tidyverse set of R packages for manipulating and plotting R data, sounds very useful. I've already used ggplot2, but the others sound great too...

Inserting pdf images in a Word doc by converting to png/jpg

I wanted to insert some pdf-format images in a Word document, but they always look fuzzy when I do this.

My colleague James suggested I convert the pdf files to a lossless image format like tiff or jpeg or png.

I was able to do this on Linux, using the ImageMagick command:
% convert -density 250 figure.pdf -quality 100 figure.png
where figure.pdf was my input pdf file, and figure.png was the output png file.

If you need to have a certain colourspace (e.g. RGB) and resolution (e.g. 300 dpi), you can use:
% convert -colorspace RGB -units PixelsPerInch -density 300 figure.pdf -quality 100 figure.png
where -density 300 specifies 300 dpi, PixelsPerInch means dpi.
Or for a jpg:
%  convert -colorspace RGB -units PixelsPerInch -density 300 figure.pdf -quality 100 figure.jpg

To check the colourspace and dpi of an image you can type:
% identify -verbose figure.jpg | more
You'll see something like this, showing it is 300 dpi and in RGB colourspace:
Image: Supplementary_Figure_1_wasExt1_5May2018.jpg
  Format: JPEG (Joint Photographic Experts Group JFIF format)
  Class: DirectClass
  Geometry: 2550x2305+0+0
  Resolution: 300x300
  Print size: 8.5x7.68333
  Units: PixelsPerInch
  Type: TrueColor
  Endianess: Undefined
  Colorspace: RGB

Problem: does your jpg look grainy?
Note that when you do this, a pdf file that looks high resolution is converted to a jpg file that looks a bit grainy. In this case it might be a helpful to convert the pdf first to a tif or gif (e.g. by doing 'Save as' in Acrobat reader and choosing 'tif' or 'gif', or by converting the pdf to a tif/gif using ImageMagick convert command), and then converting the tif/gif to a jpg using the ImageMagick convert command.

Thanks to my Dad and colleague James Cotton for helpful advice!

Information content of a GO term

I want to calculate the 'information content' (IC) of a GO term. I have a vague idea that this will tell me something about how information-rich a particular GO term is, compared to other GO terms... But how is it defined exactly?

Definition of information content of a GO term
By looking at the documentation of the GOSemSim R package, I found out the information content (IC) of a GO term is defined as the negative log probability of the term occurring in GO corpus.

The frequency of a term t is defined as: p(t)=nt′N|t′∈{t,childrenoft}

where nt′is the number of annotations with term t′, and
N is the total number of annotations in the GO corpus.
Thus the information content is defined as: IC(t)=−log(p(t)).
Here 'children of t' are all the descendants of t (see Mistry & Pavlidis).
This IC is calculated separately for terms in the 'biological process', 'molecular function' and 'cellular component' ontologies  (see Mistry & Pavlidis).

Another way of expressing this is: p(t_i) = freq(t_i)/freq(root), where 'root' is a term at the root of an ontology (for  'biological process', 'molecular function' and 'cellular component' ontologies), and  
freq(t_i) is given by:

freq(ti)=∣∣annot(ti)∣∣+∑c∈children(ti)∣∣annot(c)∣∣

(see Mistry & Pavlidis).

This means a rarely used term contains a greater amount of information. 
Mistry & Pavlidis say: 'The information content (IC) of a term is related to how often the term is applied to genes in the database, such that rarely used terms are ascribed higher IC. The IC for GO terms is monotonically decreasing as one follows the graph from a leaf terms towards the root term. Intuitively, terms low in the hierarchy are "more detailed" and impart more information about function than high-level terms such as "metabolism".'

Calculating the information content of a GO term using Python
Next I wanted to write a Python script to calculate the information content of GO terms. 

First I  made an input file with the number of genes that each GO term is assigned to in my annotation file for my species of interest. It looks something like this:
GO:0032436 1
GO:0010608 7
GO:0050577 1
GO:0005319 3
GO:0098542 1
... 
This is for GO terms in the three ontologies (biological process, molecular function, cellular component).
(Note to self: made using /nfs/helminths02/analysis/50HGP/00ANALYSES/final_GO_terms/make_files_of_go_counts.pl). 

I've written a Python script calc_information_content_for_GO_terms.py to calculate the information content for GO terms. It works like this:
% python3 calc_information_content_for_GO_terms.py go-basic.obo caenorhabditis_elegans_GO_cnts.txt caenorhabditis_elegans_IC
where go-basic.obo is your input obo (ontology hierarchy) file,
caenorhabditis_elegans_GO_cnts.txt is the file of counts of GO annotations for each GO term in yoru species of interest,
caenorhabditis_elegans_IC is the output file with information content for each GO term.

(Note to self: I find a copy of the GO ontology used by Bhavana for 50HG here: /warehouse/pathogen_wh01/users/bh4/50HGI_FuncAnnotation/go-basic.obo).

Example results of my script
The highest information content value that I calculated was 9.20 for GO:2001272, 'positive regulation of cysteine-type endopeptidase activity involved in execution phase of apoptosis'. I found in QuickGO that this is far down the GO hierarchy, as expected if it has high information content:

The lowest information content value I calculated was for GO:0005488, of 0.82. This GO term is 'binding', and is quite near the top of the GO hierarchy, as you'd expect for a low information content:

Learning R for beginners

I was just talking to my colleagues about resources for learning R for beginners, here are a couple of resources we discussed:
- 'Introductory R' by Robert Knell was highly recommended
- I like 'Kickstarting R' by Jim Lemon

Accessing EBI data using bioservices

My colleague Guillaume Salle has just pointed out that it's possible to access EBI data such as UniProt and ChEMBL using a Python module called bioservices, how cool! Here is a tutorial.

Thanks Guillaume!

Finding all occurrences of a subsequence in a sequence

I wanted to find all the occurrences of a subsequence in a genome assembly. To do this, I first tried using BLAT but it didn't find them for me (not sure why).

So I instead wrote a little Python function to print out all the positions of a subsequence in a sequence:

#====================================================================#

# find the positions of a subsequence in a sequence:

def find_positions_of_subsequence(seq, subsequence, seqname):

    still_searching = True
    start = 0
    end = len(seq) - 1
    while (still_searching == True):
        position = seq.find(subsequence, start, end)
        if position == -1:
            still_searching = False
        else:
            actual_position = position + 1
            format_string = "Found at %d in %s" % (actual_position, seqname)
            print(format_string)
            start = position + 1

    return

#====================================================================#

Python saves the day!

Trimming adapter sequences

I wanted to trim adapter sequences, specifically the sequence 'GTTTTAGGTC'

Attempt 1: Trimmomatric
I first tried Trimmomatic, with a fasta file with the forward and reverse complement of my adapter sequences, and the 'ILLUMINACLIP:/lustre/scratch118/infgen/team133/alc/000_FUGI_PatrickCRISPR/OmegaGenewhiz/adaptor.fa:2:40:0' option. However, I found that this just removes the sequences with adapter, it doesn't trim them. This is also described here.

Attempt 2: Cutadapt
I then tried cutadapt. With the following fasta file:
>seq1
GTTTTAGGTCGTTATCGTGTA
>seq2
TACACGATAACGACCTAAAAC 

I was able to trim adapters using:
% cutadapt -g GTTTTAGGTC my.fasta -a GACCTAAAAC
where -g cuts sequences off the 5' end, and -a off the 3' end. 
It cuts off adapter with at most 10% errors in single-end read mode.

This gave me:
>seq1
GTTATCGTGTA
>seq2
TACACGATAAC
Hurray!

Something strange I noticed about cutadapt:
I have been using cutadapt to trim some adaptors, and noticed that it totally removed some of my sequences.

For example, with this fastq input: (in tmp.fastq)
@M03558:259:000000000-BH588:1:1101:4913:6565 2:N:0:TTTGTA
ACTGACCCTCAGCAATCTTAAACTTCTTAGACGAATCACCAGAACGGAAAACATCCTTCATAGAAATTTCACGCGGCGGCAAGTTGCCATACAAAACAGGGTCGCCAGCAATATCGGTATAAGTCAAAGCACCTTTAGCGTTAAGGTACTGAATCTCTTTAGTCGCAGTAGGCGGAAAACGAACAAGCGCAAGAGTAAACATAGTGCCATGCTCAGGAACAAAGAAACGCGGCACAGAATGTTTATAGGTC
+
AAAAAFFFFCFCGGGGGGGGGGHHHHHHGFHHGEEGEHHFHHFFHGEGFGHFHHHHHHHHHHGHFHHHHHHHEHGCFGGGGGGHHHHHHHHHHHHHGHGHGHHGGGGGGHHGFHEGGCEFFHHGGGHHHEFGHEHHHGGGAE?EHHGDHHFGHHHHHEHHHH:C.A@DCGHBGA-?BGGGCGEGCFFFFFFFFEF;FFFEF/B/:BFEF;/:;BB;BFFFFFF/BFFFAAFA-;-B.FFBF;/BFBF/9/:

I then ran cutadapt like this:
% cutadapt -g GTTTTAGGTC -a GACCTAAAAC tmp.fastq
and it removed the whole input sequence.
I can't see any matches to the adaptors in that input though. When I ran blastn between that read and the NCBI database, the sequence was found to be PhiX174. Is that why cutadapt removed it, or was there some other reason?

Update: 7-Dec-2017: I sent an email about this to the cutadapt developer, Marcel Martin, and received a very helpful reply:

In short: The 5' adapter is found at the very end (3') of the read with one 

error.

From the statistics report that cutadapt prints, you can see that it finds 

the 5' adapter with one error. And if you run the above command without 

specifying the 3' adapter with -a GACCTAAAAC, you get the same result, so 

it must be the 5' adapter. When cutadapt finds a 5' adapter, it removes 

the adapter sequnce itself and everything preceding it. And in your case, 

the read ends with 'GTTTATAGGTC', which is the 5' adapter sequence with an 

extra 'A' inserted between one of the 'T' nucleotides.
Although cutadapt doesn’t print it, the alignment likely looks like this:

...GTTTATAGGTC (read end)
   GTTT-TAGGTC (adapter)

This is possibly not the result you were expecting, but cutadapt is doing 

what it is supposed to.

You have a few options:
- Accept this result
- If you think allowing insertions is the problem, use --no-indels to
disallow insertions and deletions.
- If you want to prevent the 5' adapter from occurring within the read at 

all, you can specify the adapter as -g XXGTTTTAGGTC. The 'X' is always 

counted as a mismatch, so if you have more X than allowed mismatches, 

the adapter is forced to occur at the 5' end (overlapping it).