On Jul 17, 2008, at 9:47 AM, Krys Kelly wrote:
I have inherited a pipeline for Solexa sequence data ... This
seems pretty
crude
could anyone point me in the direction of some code?
[ This may now be silly, but I will send it anyway. ]
On (quasi) UNIX, you can use something like the script below, which
elides the
largest initial portion of our 3' adapter and then the largest
terminal portion
of our 5' adapter, minimal length 3 in both cases, from a data file
where the
read (possibly containing ".") was the 5th field.
To allow 1 error, (using awk, for example) you can blow each pattern
in there
up into a regular expression as in this sed edit command (the \'s
don't belong
but were necessary for some reason):
's,.CGTATGCCGTCTTCTGCTTG$\|T.GTATGCCGTCTTCTGCTTG$\|
TC.TATGCCGTCTTCTGCTTG$\|TCG.ATGCCGTCTTCTGCTTG$\|TCGT.TGCCGTCTTCTGCTTG$
\|TCGTA.GCCGTCTTCTGCTTG$\|TCGTAT.CCGTCTTCTGCTTG$\|
TCGTATG.CGTCTTCTGCTTG$\|TCGTATGC.GTCTTCTGCTTG$\|TCGTATGCC.TCTTCTGCTTG$
\|TCGTATGCCG.CTTCTGCTTG$\|TCGTATGCCGT.TTCTGCTTG$\|
TCGTATGCCGTC.TCTGCTTG$\|TCGTATGCCGTCT.CTGCTTG$\|TCGTATGCCGTCTT.TGCTTG$
\|TCGTATGCCGTCTTC.GCTTG$\|TCGTATGCCGTCTTCT.CTTG$\|
TCGTATGCCGTCTTCTG.TTG$\|TCGTATGCCGTCTTCTGC.TG$\|TCGTATGCCGTCTTCTGCT.G$
\|TCGTATGCCGTCTTCTGCTT.$,>,'
In our data, the whole 3' adapter may occur at the 5' end, making the
read
worthless, or in the interior!
All of this probably translates into R, but this was quicker for me.
-Harris
--------------
#!/bin/sh
DIR=Processed
fgrep -v . |
awk '{print $5}' |
sed \
-e 's,TCGTATGCCGTCTTCTGCTTG$,>,' \
-e 's,TCGTATGCCGTCTTCTGCTT$,>,' \
-e 's,TCGTATGCCGTCTTCTGCT$,>,' \
-e 's,TCGTATGCCGTCTTCTGC$,>,' \
-e 's,TCGTATGCCGTCTTCTG$,>,' \
-e 's,TCGTATGCCGTCTTCT$,>,' \
-e 's,TCGTATGCCGTCTTC$,>,' \
-e 's,TCGTATGCCGTCTT$,>,' \
-e 's,TCGTATGCCGTCT$,>,' \
-e 's,TCGTATGCCGTC$,>,' \
-e 's,TCGTATGCCGT$,>,' \
-e 's,TCGTATGCCG$,>,' \
-e 's,TCGTATGCC$,>,' \
-e 's,TCGTATGC$,>,' \
-e 's,TCGTATG$,>,' \
-e 's,TCGTAT$,>,' \
-e 's,TCGTA$,>,' \
-e 's,TCGT$,>,' \
-e 's,TCG$,>,' \
-e 's,^GTTCAGAGTTCTACAGTCCGACGATC,<,' \
-e 's,^TTCAGAGTTCTACAGTCCGACGATC,<,' \
-e 's,^TCAGAGTTCTACAGTCCGACGATC,<,' \
-e 's,^CAGAGTTCTACAGTCCGACGATC,<,' \
-e 's,^AGAGTTCTACAGTCCGACGATC,<,' \
-e 's,^GAGTTCTACAGTCCGACGATC,<,' \
-e 's,^AGTTCTACAGTCCGACGATC,<,' \
-e 's,^GTTCTACAGTCCGACGATC,<,' \
-e 's,^TTCTACAGTCCGACGATC,<,' \
-e 's,^TCTACAGTCCGACGATC,<,' \
-e 's,^CTACAGTCCGACGATC,<,' \
-e 's,^TACAGTCCGACGATC,<,' \
-e 's,^ACAGTCCGACGATC,<,' \
-e 's,^CAGTCCGACGATC,<,' \
-e 's,^AGTCCGACGATC,<,' \
-e 's,^GTCCGACGATC,<,' \
-e 's,^TCCGACGATC,<,' \
-e 's,^CCGACGATC,<,' \
-e 's,^CGACGATC,<,' \
-e 's,^GACGATC,<,' \
-e 's,^ACGATC,<,' \
-e 's,^CGATC,<,' \
-e 's,^GATC,<,' \
-e 's,^ATC,<,' | tr -d '<>' |
awk 'length>15' | sort | uniq -c |
awk '{print >DIR "/" length($2)}' DIR=$DIR
_______________________________________________
Bioc-sig-sequencing mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing
