new script gmapRescue.sh

bin/gmapRescue.sh

+#!/bin/bash
+##SBATCH --array=0-7
+##SBATCH -n1
+##SBATCH --mail-user=helene.rimbert@inrae.fr
+##SBATCH --mail-type=ALL
+
+# for use with slurm
+#chroms=($(echo TraesCS{1..7} TraesCSU))
+#chrom=${chroms[$SLURM_ARRAY_TASK_ID]}
+
+# get the list of query and target file for gmap
+
+inputDir=$1
+blastdb=$2
+
+querylist=($(find  $inputDir  -type f -name 'query.fasta'))
+targetlist=($(find $inputDir -type f -name 'target.fasta'))
+
+numquery=${#querylist[@]}
+numtarget=${#targetlist[@]}
+#blastdb='/banks/blast/IWGSCv1.1_repr_mrna'
+
+echo "STARTING AT `date` FOR DATASET $inputDir"
+echo "found $numquery query fasta file"
+echo "found $numtarget target fasta file"
+
+
+for i in "${querylist[@]}";
+do
+	rootdir=`dirname $i`
+	# extract the name of the gene
+        geneid=$(cat $rootdir'/query.fasta'|fgrep '>'|perl -ne 'print $1 if /^\>(TraesCS.+)\ coords/')
+	echo "\n\n================================================"
+        echo " gene is $geneid"
+	echo " datadir is $rootdir"
+
+	echo "1. creating gmap db for dataset $rootdir"
+
+	if [ ! -d "$rootdir/query" ]
+	then
+		gmap_build -D $rootdir -d query $rootdir'/query.fasta' &
+		pid=$!
+		gmap_build -D $rootdir -d target $rootdir'/target.fasta' &
+		pid2=$!
+		wait $pid $pid2
+	fi
+
+	echo "2. extract mrna"
+	mrnaid=$geneid'.1'
+
+	if [ ! -s "$rootdir/query.mrna.fasta" ]
+	then
+		blastdbcmd -db $blastdb -entry $mrnaid > $rootdir'/query.mrna.fasta'
+		
+	fi
+
+	echo "3. run gmap on target and query"
+	if [ ! -s "$rootdir/gmap.target.gff3" ]
+	then 
+		gmap -d target -D $rootdir -f 2 $rootdir'/query.mrna.fasta' > $rootdir'/gmap.target.gff3' &
+		pid=$!
+		gmap -d query -D $rootdir -f 2 $rootdir'/query.mrna.fasta' > $rootdir'/gmap.query.gff3' &
+		pid2=$!
+		blastn -query $rootdir'/query.fasta' -subject $rootdir'/target.fasta' -outfmt 6 -out $rootdir'/blast_m8.tab' &
+		pid3=$!
+		wait $pid $pid2 $pid3
+	fi
+	rm -rf $rootdir'/query' $rootdir'/target'
+
+	echo "4. indexing the query and target genomic"
+	if [ ! -s " $rootdir/query.fasta.fai" ]
+	then
+		samtools faidx $rootdir'/query.fasta'
+		samtools faidx $rootdir'/target.fasta'
+	else
+		echo " index already done"
+	fi
+
+	echo "5. running gffread to extract the pep sequence"
+
+	if [ ! -s " $rootdir/target.pep.faa" ]
+	then
+		gffread -x $rootdir'/target.cds.fasta' -y $rootdir'/target.pep.faa' -g $rootdir'/target.fasta' $rootdir'/gmap.target.gff3'
+		gffread -x $rootdir'/query.cds.fasta' -y $rootdir'/query.pep.faa' -g $rootdir'/query.fasta' $rootdir'/gmap.query.gff3'
+	else
+		echo "pep extraction already done"
+	fi
+
+	echo "6. running diff to check if the pep sequences are different or not"
+	diffresult=$(diff $rootdir'/target.pep.faa' $rootdir'/query.pep.faa')
+
+	if [ ! -s "$rootdir/blat.txt" ]
+	then
+		echo "RESULTS no blat results here: no extraction to make"
+		touch $rootdir/UNMAPPED
+	else
+		if [ "$diffresult" != "" ]
+		then
+			echo "RESULTS pep of gene $geneid has changed"
+			echo $diffresults > $rootdir/CDS_CHANGED
+
+		else
+			echo "RESULTS pep of gene $geneid stays the same"
+			touch $rootdir/CDS_OK
+		fi
+	fi
+
+done
+
+echo "ENDED AT `date`"