Help with mapping sgRNA reads from CRISPR/Cas9 screen

[N00bSeq]

Hi!

I am somewhat of a beginner and will soon receive data from a CRISPR-screen, and would like to get familiar with the workflow. For now, I am trying to map and normalize reads from public data.

The reads I am using are obtained from here:

http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM1715941

The reference sgRNA sequences corresponding to this study (as far as I can tell) (GeCKOv1 library) were obtained here:

https://sourceforge.net/p/mageck/wiki/libraries/

I have made a fasta file from the GeCKOv1.txt file and indexed it with bowtie 1.0.0 in the following way:

awk '{print ">"$1"\n"$2}' GeCKOv1.txt > GeCKOv1.fasta
bowtie-build GeCKOv1.fasta GeCKOv1

However, when attempting to map the reads from the sample SRR2071309.fastq to this reference (as follows) 0 reads are mapped:

bowtie -t GeCKOv1 SRR2071309.fastq SRR2071309.GeCKOv1.map

I suspect that this may be because the reads need to be trimmed somehow, since sgRNAs are generally only about 20 bp and the reads are 51 bp. Therefore, I printed the first few lines for the fastq-file to identify any recurring sequences in the beginning or end:

@SRR2071309.1 DH1DQQN1:414:H9PTDADXX:1:1101:1913:2190 length=51
ACTGATTCTTGTGGAAAGGACGAAACACCGCGTCGGGATGCACCAGCTCCG
+SRR2071309.1 DH1DQQN1:414:H9PTDADXX:1:1101:1913:2190 length=51
DDDDDDDIIIIEFFIEIII>EEIIIIIIDIIDIIDDDDDDDDDAA@@AAAA
@SRR2071309.2 DH1DQQN1:414:H9PTDADXX:1:1101:1832:2196 length=51
ACTGTTCTTGTGGAAAGGACGAAACACCGTAGATACCCAGAACGCCTTCGT
+SRR2071309.2 DH1DQQN1:414:H9PTDADXX:1:1101:1832:2196 length=51
FFHHCFHJJJIIJJJJJJJJJJJJJJJJJHIIIJJJJJJJJJJJJHHHFFF
@SRR2071309.3 DH1DQQN1:414:H9PTDADXX:1:1101:2186:2239 length=51
ACTGATTCTTGTGGAAAGGACGAAACACCGCACACAGAAGAGCTCCTGGCG

I made a guess that "ACTGATTCTTGTGGAAAGGACGAAACACCG" would be the sequence that should be trimmed. Therefore I tried cutadapt:

cutadapt -a ACTGATTCTTGTGGAAAGGACGAAACACCG -o SRR2071309.trimmed.fastq SRR2071309.fastq

However, the output file (around 98% reads trimmed according to the log) looked like this:

@SRR2071309.1 DH1DQQN1:414:H9PTDADXX:1:1101:1913:2190 length=51

+SRR2071309.1 DH1DQQN1:414:H9PTDADXX:1:1101:1913:2190 length=51

@SRR2071309.2 DH1DQQN1:414:H9PTDADXX:1:1101:1832:2196 length=51

+SRR2071309.2 DH1DQQN1:414:H9PTDADXX:1:1101:1832:2196 length=51

@SRR2071309.3 DH1DQQN1:414:H9PTDADXX:1:1101:2186:2239 length=51

I am completely lost as to how I should proceed to be able to map these reads to the reference. Is anyone familiar with this?

[cmbetts]

I don't use cutadapt, so I'm not familiar with it's options, but my guess would be that it's trimming everything after your adapter sequence, leaving nothing behind.
In most NGS applications, the need for adapter trimming is when the library inserts are shorter than the read lengths, so the read goes through to the adapter on the opposite side. In those cases, everything after the adapter is garbage sequence that needs to be trimmed before further analysis. You data has a funky structure where the adapter is the 5' most end of the reads, so when the program runs expecting typical library structure it discards just about everything.
You need to do the trimming in a way that's aware the adapter is at the start of the read and that you want to keep the trailing sequence. I'm not sure which programs have this option, I've been using the BBMAP tools lately, and they generally have options galore. Another possible option would be to reverse compliment your reads so that the structure looked like a typical library, and then just about any trimming program could handle it no problem.

[N00bSeq]

You are right. That was an issue. I instead tried trimming a fix number of bases from the 5' end with another tool

seqtk trimfq -b 30 SRR2071309.fastq > SRR2071309.setk.trimmed.fastq

Then mapping again:

bowtie -t GeCKOv1 SRR2071309.setk.trimmed.fastq SRR2071309.GeCKOv1.setk.trimmed.map

However, this also produced 0 alignments.

Cannot figure this out.

[N00bSeq]

It must be something more fundamental that I'm doing wrong here. Now I also tried to make use of the info on https://sourceforge.net/p/mageck/wiki/Home/, up to the part about trimming:

I downloaded ERR376998, trimmed the 23 first bases according to the instructions and mapped them to the "yusa_library" reference, again with 0 % mapped reads. Does the fasta reference (converted from the yusa_library.csv file from the above link) seem to look OK?

>chr9:21066101-21066124
CCTTCTCAATACCCATTCG
>chr9:21066102-21066125
CTTCTCAATACCCATTCGC
>chr9:21066273-21066296
GACCTGGGTGGCCGGCGAA
>chr9:21066274-21066297
ACCTGGGTGGCCGGCGAAG
>chr9:21066926-21066949
CGCGCCACGCACCTGATAT

Or what else could be the problem?

[Dario1984]

The U6 promoter sequence is not a fixed number of bases from the 5' end. In my dataset, I see the start of the promoter frequently between positions 10 and 18.

A better method is to find the start and end positions of the U6 promoter sequence in your reads and also the plasmid sequence (i.e. GTTTT...) and extract the sequence in between. Once you have extracted the sequence in between you need to filter the sequences to be between 19 and 21 bases long. Some of the extracted sequences can be 30 or 40 nucleotides long, because two different gRNAs sometimes dimerise. For example, a 75 base read from my dataset is:

TAACAATGGTCTTGTGGAAAGGACGAAACACCA[COLOR="Blue"]AGGTAACA[/COLOR]A[COLOR="Red"]ATGCGTGCGAGCCG[/COLOR]GTTTTAGAGCTAGAAATAG
                             TTCT[COLOR="Blue"]AGGTAACA[/COLOR]GGAGTTGT (TMTC3)
                                    TTCCTC[COLOR="Red"]ATGCGTGCGAGCCG[/COLOR] (SIVA1)

See how the read contains (parts of) two gRNAs ? Also the U6 promoter ends in A, not G, for this read. I typically see 80% ending in CACCG and 17% ending in CACCA. This creates additional complexity for finding the end of the U6 promoter.

There's a lot of variability in CRISPR experiments. You need to find the preprocessing parameters suited to your dataset.