/usr/share/doc/last-align/examples/last-bisulfite-paired.sh is in last-align 393-1.
This file is owned by root:root, with mode 0o755.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 | #! /bin/sh
# Align paired bisulfite-converted DNA reads to a genome.
# This assumes that reads1.fastq are all from the converted strand
# (i.e. they have C->T conversions) and reads2.fastq are all from the
# reverse-complement (i.e. they have G->A conversions).
# "GNU parallel" needs to be installed.
[ $# -eq 4 ] || {
cat <<EOF
Typical usage:
lastdb -w 2 -u bisulfite_f.seed my_f mygenome.fa
lastdb -w 2 -u bisulfite_r.seed my_r mygenome.fa
$(basename $0) my_f my_r reads1.fastq reads2.fastq > results.maf
EOF
exit 2
}
# Try to get the LAST programs into the PATH, if they aren't already:
PATH=$PATH:$(dirname $0)/../src:$(dirname $0)/../scripts
tmp=${TMPDIR-/tmp}/$$
trap 'rm -f $tmp.*' EXIT
cat > $tmp.fmat << 'EOF'
A C G T
A 6 -18 -18 -18
C -18 6 -18 3
G -18 -18 6 -18
T -18 -18 -18 3
EOF
cat > $tmp.rmat << 'EOF'
A C G T
A 3 -18 -18 -18
C -18 6 -18 -18
G 3 -18 6 -18
T -18 -18 -18 6
EOF
cat > $tmp.script << 'EOF'
t=$1.$$
lastal -p $1.fmat -s1 -Q1 -e120 -i1 "$2" "$4" > $t.t1f
lastal -p $1.rmat -s0 -Q1 -e120 -i1 "$3" "$4" > $t.t1r
last-merge-batches.py $t.t1f $t.t1r > $t.t1
rm $t.t1f $t.t1r
lastal -p $1.fmat -s0 -Q1 -e120 -i1 "$2" "$5" > $t.t2f
lastal -p $1.rmat -s1 -Q1 -e120 -i1 "$3" "$5" > $t.t2r
last-merge-batches.py $t.t2f $t.t2r > $t.t2
rm $t.t2f $t.t2r
last-pair-probs.py -m0.1 $t.t1 $t.t2 |
perl -F'(\s+)' -ane '$F[12] =~ y/ta/CG/ if /^s/ and $s++ % 2; print @F'
rm $t.t1 $t.t2
EOF
# Convert C to t, and all other letters to uppercase:
perl -pe 'y/Cca-z/ttA-Z/ if $. % 4 == 2' "$3" | split -l400000 -a5 - $tmp.1
# Convert G to a, and all other letters to uppercase:
perl -pe 'y/Gga-z/aaA-Z/ if $. % 4 == 2' "$4" | split -l400000 -a5 - $tmp.2
parallel --gnu --xapply sh $tmp.script $tmp "$1" "$2" ::: $tmp.1* ::: $tmp.2*
|