SEQanswers

Go Back   SEQanswers > Applications Forums > RNA Sequencing



Similar Threads
Thread Thread Starter Forum Replies Last Post
CRISPR/Cas9 sequencing muzz56 Illumina/Solexa 1 03-17-2015 12:01 PM

Reply
 
Thread Tools
Old 03-09-2016, 07:23 AM   #1
N00bSeq
Member
 
Location: Sweden

Join Date: Mar 2014
Posts: 12
Default Help with mapping sgRNA reads from CRISPR/Cas9 screen

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/quer...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:

Code:
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:

Code:
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:

Code:
@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:

Code:
cutadapt -a ACTGATTCTTGTGGAAAGGACGAAACACCG -o SRR2071309.trimmed.fastq SRR2071309.fastq
However, the output file (around 98% reads trimmed according to the log) looked like this:

Code:
@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?
N00bSeq is offline   Reply With Quote
Old 03-09-2016, 09:17 AM   #2
cmbetts
Senior Member
 
Location: Bay Area

Join Date: Jun 2012
Posts: 114
Default

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.
cmbetts is offline   Reply With Quote
Old 03-09-2016, 09:34 AM   #3
N00bSeq
Member
 
Location: Sweden

Join Date: Mar 2014
Posts: 12
Default

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

Code:
seqtk trimfq -b 30 SRR2071309.fastq > SRR2071309.setk.trimmed.fastq
Then mapping again:

Code:
bowtie -t GeCKOv1 SRR2071309.setk.trimmed.fastq SRR2071309.GeCKOv1.setk.trimmed.map
However, this also produced 0 alignments.

Cannot figure this out.
N00bSeq is offline   Reply With Quote
Old 03-09-2016, 09:54 AM   #4
N00bSeq
Member
 
Location: Sweden

Join Date: Mar 2014
Posts: 12
Default

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?

Code:
>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?
N00bSeq is offline   Reply With Quote
Old 12-13-2016, 09:05 PM   #5
Dario1984
Senior Member
 
Location: Sydney, Australia

Join Date: Jun 2011
Posts: 166
Default

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:

Code:
TAACAATGGTCTTGTGGAAAGGACGAAACACCAAGGTAACAAATGCGTGCGAGCCGGTTTTAGAGCTAGAAATAG
                             TTCTAGGTAACAGGAGTTGT (TMTC3)
                                    TTCCTCATGCGTGCGAGCCG (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.
Dario1984 is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 02:24 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO