SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Extracting sub-sequences from a larger DNA sequence using a file of coordinates gwilymh Bioinformatics 3 03-04-2015 02:54 AM
Extract sequences from a FASTQ file based on another file caputcastellae Bioinformatics 3 08-14-2014 01:39 PM
Fastq file format for paired end sequences rozitaa Bioinformatics 9 07-03-2013 03:25 AM
Sorting fastq by primers, then searching by sequence (with mismatches) jme Bioinformatics 0 01-18-2012 09:25 AM
extracting reads from a large fastq file Wallysb01 Bioinformatics 23 08-08-2011 01:43 PM

Reply
 
Thread Tools
Old 10-26-2015, 12:29 PM   #1
rtleenay
Junior Member
 
Location: North Carolina

Join Date: Oct 2015
Posts: 4
Default Searching for and extracting sequences within a fastQ file

Hello,

So I'm trying to extract information from a nucleotide library I have sequenced on an Illumina MiSeq.

In short, I amplified an ~150bp region, which contains a 4 nucleotide library about 70 nucleotides into the read. This means that about 90% of each read is unnecessary, and all I need is the small 4 bp region.

Previously, I have grabbed out my region of interest using the following code:

grep 'TTCATTAAAAATTGAATTGACATTAACCTATAAAAATAGGCGTCGAGGCCCTTTCGTCTTC[TCAG][TCAG][TCAG][TCAG]GTCGAGTGCA' sample.fastq | cut -c 62-65 | sort | uniq -c | sort -nr | less > sample.txt

This code cuts out the library region of interest based on the character number (see cut -c) and then sorts and counts how many times that sequence occurred.

My problem is that now I have done this again, except with a slightly modified sequence. This time, there are a few homopolmeric regions on either side of my library, so occasionally they are miscalled as one nucleotide longer or shorter.

The sequence looks like this: GGGNNNNGGG

Instead of cutting by character, I would just like to cut out this region from each read, and then put it through the rest of my counting pipeline. However, I cannot seem to find a good way to cut out only this sequence from within the fastq file.

Does anyone have any suggestions to fix this?

Thank you!

-rtleenay
rtleenay is offline   Reply With Quote
Old 10-26-2015, 01:23 PM   #2
mastal
Senior Member
 
Location: uk

Join Date: Mar 2009
Posts: 667
Default

changing to something like:

grep -P 'GGG{min, max}[ACGT]{4}GGG{min,max} sample.fastq |

might work.

where min and max could be set to one less (for min), or one more (max) than the number of nucleotides in your homopolymer stretch.
mastal is offline   Reply With Quote
Old 10-27-2015, 05:31 AM   #3
rtleenay
Junior Member
 
Location: North Carolina

Join Date: Oct 2015
Posts: 4
Default

But then once I do use grep to grab that line, what command (or tool) can I use to grab that stretch by its sequence instead of by the character number?
rtleenay is offline   Reply With Quote
Old 10-27-2015, 06:45 AM   #4
mastal
Senior Member
 
Location: uk

Join Date: Mar 2009
Posts: 667
Default

In perl, you could capture the sequence that matches the [ACGT]{4} using regular expressions, and print it to a file for downstream processing.
mastal 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 09:43 AM.


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