Skip to content
Snippets Groups Projects
Commit a4418dc5 authored by Helene Rimbert's avatar Helene Rimbert
Browse files

remove useless scripts extractpepandcheck.sh recalcGmapRescue.sh

parent d553b5f7
No related branches found
No related tags found
No related merge requests found
#!/bin/bash
#SBATCH --array=0-7
#SBATCH -n1
#SBATCH --mail-user=helene.rimbert@inrae.fr
#SBATCH --mail-type=ALL
#SBATCH --export=ALL
chroms=($(echo TraesCS{1..7} TraesCSU))
chrom=${chroms[$SLURM_ARRAY_TASK_ID]}
# get the list of query and target file for gmap
querylist=($(find 3_mapping_$chrom* -type f -name 'query.fasta'))
targetlist=($(find 3_mapping_$chrom* -type f -name 'target.fasta'))
numquery=${#querylist[@]}
numtarget=${#targetlist[@]}
echo "STARTING AT `date` FOR CHROM $chrom"
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. 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 "2. 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 "3. running diff to check if the pep sequences are different or not"
diffresult=$(diff $rootdir'/target.pep.faa' $rootdir'/query.pep.faa')
if [ "$diffresult" != "" ]
then
echo "RESULTS pep of gene $geneid has changed"
echo $diffresults > $rootdir/CDS_CHANGED
elif [ ! -s "$rootdir/blat.txt" ]
then
echo "RESULTS no blat results here: no extraction to make"
touch $rootdir/UNMAPPED
else
echo "RESULTS pep of gene $geneid stays the same"
touch $rootdir/CDS_OK
fi
done
echo "ENDED AT `date`"
#!/bin/bash
##SBATCH --array=0-7
##SBATCH -n1
##SBATCH --mail-user=helene.rimbert@inrae.fr
##SBATCH --mail-type=ALL
inputDir=$1
inputDescription=$2
finalGFF=$3
bindir=$4
export TRANSFERTANNOT_PATH=$bindir
querylist=($(find $inputDir -type f -name 'gmap.target.gff3'))
numquery=${#querylist[@]}
echo "STARTING AT `date` FOR DATASET $inputDir"
echo "found $numquery query file"
for inputgff in "${querylist[@]}";
do
rootdir=`dirname $inputgff`
# 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"
# extract the target chrom, start and stop
chrom=$(fgrep '>' $rootdir'/target.fasta'|sed "s/>//"|cut -d ' ' -f1|cut -d '_' -f2)
start=$(fgrep '>' $rootdir'/target.fasta'|sed "s/>//"|cut -d ' ' -f1|cut -d '_' -f3)
stop=$(fgrep '>' $rootdir'/target.fasta'|sed "s/>//"|cut -d ' ' -f1|cut -d '_' -f4)
echo " Target coords: $chrom from $start - $stop"
# check the status of the CDS: OK or CHANGED or UNKNOWN
mappingStatus='UNMAPPED'
cdsStatus='CDS_OK'
if [ -f "$rootdir/CDS_OK" ]
then
echo " CDS has not changed between target and reference"
cdsStatus='CDS_OK'
elif [ -f "$rootdir/CDS_CHANGED" ]
then
echo " WARNING : CDS has changed between target and reference"
cdsStatus='CDS_CHANGED'
elif [ -f "$rootdir/CDS_UNKNOWN" ]
then
echo " WARNING : we do not know if the CDS has changed or not"
cdsStatus='CDS_UNKNOWN'
else
echo " Assume gene is UNKNOWN"
echo " need to run gmap on entire genome"
cdsStatus='CDS_UNKNOWN'
continue
fi
echo " 1. convert gmap coordinates on new reference"
eval $TRANSFERTANNOT_PATH'/bin/convertGmapCoords.py' -G $inputgff -s $start -e $stop -c $cdsStatus -o $rootdir'/recalc.gff3' -r $chrom -g $geneid -m 'GMAP_ON_TARGET' &> $rootdir'/recalc.log'
echo " 2. sort gff"
eval gt gff3 -sort -retainids -tidy $rootdir'/recalc.gff3' 1> $rootdir'/recalc_clean.gff3'
done
# merge All recalc_clean.gff3 files
gfflist=($(find $inputDir -type f -name 'recalc_clean.gff3'))
gt gff3 -sort -tidy -retainids "${gfflist[@]}" 1> $finalGFF 2> $finalGFF.gt.log
echo "ENDED AT `date`"
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment