search multiple sequences in multiple fastqs
Hi all!
I have a script which goes through a bunch of fastq files and searches for specific sequences found in a text file called test_text.txt, then outputs a fastq file with only those reads. The reads in my fastq files are all between 15-30bp long, however when I look at the output fastq files, I am getting reads of 52 bp long.
Does anyone know why this would happen? code is below:
for fastq in `ls my_files`;do
for i in `cat test_text.txt`; do zgrep -B 1 -A 2 $i $fastq >>filtered_fastq/${fastq}.filtered_reads.fastq; done;
done
Thanks for any input.
0 answers
No answers yet.
Log in to answer this question.
More posts like this
-
Biostrings extract read
written by marco.barr •
-
Help with the reads filtering process
written by valentinavan •
-
bioconda sabre reads shorter
written by JG.Schonthaler •
-
Size difference using gunzip-cat-gzip and zcat
written by dazhudou1122
-
Varying assembly statistics: Abyss & SPades
written by max_19
-
max length in cutadapt
written by max_19
-
STAR or Bowtie for small RNA seq?
written by max_19
-
Iterative trimming using bowtie
written by ancient_learner
-
Converting a bunch or SRA files using fastq-dump split-files
written by Mitra •
-
Merge fastq file from same lane and start bowtie alignment
written by giorgiocasaburi
What have you done so far to convince yourself/check the reads are all 15-30bp?
I can't see an obvious mistake at first glance, so are you sure there are no 52 bp reads in the input file?
I used cutadapt to set min/max length for all sequences (cutadapt -a $adapter3 -f fastq -m 15 -M 30). also when I do fastqc on the file, my sequence length in the basic stats section shows '15-30'
Make sure; do
awk '{print length}' <filename> | sort -n | tail. The last number should be the longest sequence length in your file. Remove the headers before if you think they're longer than 30bp (i.e.grep -v ">" <filename> | awk '{print length}' | sort -n | tailThank you! Tried this and got 54. This is odd! does this mean something is wrong with my cutadapt step? and why does FASTQC report 15-30bp for sequence length
cutadapt should get rid of reads > 30 in length in that case. Are you sure that you have the trimmed file in
my_files, and not the original file you had before running cutadapt?By the way, the
-fswitch in grep allows you to provide a file of patterns to search; e.g.grep -f <file_with_patterns> <file_to_search_in>, so possibly your could get rid of your inner for loop.Good to know, thank you!