Skip to content

P
pharmAnnot
Commit 0208c5db authored by Simon Verdez's avatar Simon Verdez
Browse files
Options

pipeline epilepsie

parent 99b711dc
Show whitespace changes
Inline Side-by-side
Showing with 32 additions and 308 deletions
Dup_Data.sh 0 → 100755
View file @ 0208c5db
#!/bin/bash
if [ ! -e /work/gad/sv347413/epilepsie/vcf_tmp/$i*.vcf ]
then
for i in $(cat /work/gad/sv347413/HLAscan/liste_patients_epileptiques.txt)
do
cp /work/gad/shared/vcf/individu/$i/*.raw.vcf /work/gad/sv347413/epilepsie/vcf_tmp/.
done
fi
if [ ! -e /work/gad/sv347413/epilepsie/tsv_tmp/$i*.tsv ]
then
for i in $(cat /work/gad/sv347413/HLAscan/liste_patients_epileptiques.txt)
do
cp /work/gad/shared/vcf/individu/$i/*.cnv.annot.tsv /work/gad/sv347413/epilepsie/tsv_tmp/.
done
fi
Master.sh
View file @ 0208c5db
...@@ -64,7 +64,6 @@ REF=`grep "reference" $PHARMCONFIG | cut -f2` ...@@ -64,7 +64,6 @@ REF=`grep "reference" $PHARMCONFIG | cut -f2`
PYTHONPATH=$PIPELINEBASE/common:/work/gad/shared/bin/lib/python_2.7/lib/python2.7/site-packages:/work/gad/shared/bin/miniconda2/lib/python2.7/site-packages/ PYTHONPATH=$PIPELINEBASE/common:/work/gad/shared/bin/lib/python_2.7/lib/python2.7/site-packages:/work/gad/shared/bin/miniconda2/lib/python2.7/site-packages/
export PYTHONPATH export PYTHONPATH
# Sample list # Sample list
samples=`find -L $ANALYSISDIR/ -maxdepth 1 -mindepth 1 -type d | xargs -l basename` samples=`find -L $ANALYSISDIR/ -maxdepth 1 -mindepth 1 -type d | xargs -l basename`
...@@ -75,16 +74,18 @@ then ...@@ -75,16 +74,18 @@ then
fi fi
find -L $ANALYSISDIR/ -maxdepth 1 -mindepth 1 -type d | xargs -l basename >> $ANALYSISDIR/batch.list find -L $ANALYSISDIR/ -maxdepth 1 -mindepth 1 -type d | xargs -l basename >> $ANALYSISDIR/batch.list
HLAscan:
# HLA genotyping # HLA genotyping
echo "### Execute HLAscan ###" echo "### Execute HLAscan ###"
echo "Start : $(date +"%F_%H-%M-%S")" echo "Start : $(date +"%F_%H-%M-%S")"
for currentSample in $samples for currentSample in $samples
do do
echo "Command : qsub -pe smp 2 -N HLAscan_${currentSample} -o $ANALYSISDIR/$currentSample/logs/ -e $ANALYSISDIR/$currentSample/logs/ -v INPUTFILE=$ANALYSISDIR/$currentSample/$currentSample.dedup.bam,OUTPUTFILE=$ANALYSISDIR/$currentSample/HLAscan,GeneList=$GL,LOGFILE=$ANALYSISDIR/$currentSample/logs/process_HLAscan.$logbasename.log,CONFIGFILE=$CONFIGFILE $PIPELINEBASE/haplo_scan.sh" echo "Command : qsub -pe smp 1 -N HLAscan_${currentSample} -o $ANALYSISDIR/$currentSample/logs/ -e $ANALYSISDIR/$currentSample/logs/ -v INPUTFILE=$ANALYSISDIR/$currentSample/$currentSample.dedup.bam,OUTPUTFILE=$ANALYSISDIR/$currentSample/HLAscan,GeneList=$GL,LOGFILE=$ANALYSISDIR/$currentSample/logs/process_HLAscan.$logbasename.log,CONFIGFILE=$CONFIGFILE $PIPELINEBASE/haplo_scan.sh"
qsub -pe smp 1 -N HLAscan_${currentSample} -o $ANALYSISDIR/$currentSample/logs/ -e $ANALYSISDIR/$currentSample/logs/ -v BAM=$ANALYSISDIR/$currentSample/$currentSample.dedup.bam,SORTIE=$ANALYSISDIR/$currentSample/HLAscan,GeneList=$GL,LOGFILE=$ANALYSISDIR/$currentSample/logs/process_HLAscan.$logbasename.log,CONFIGFILE=$CONFIGFILE $PIPELINEBASE/haplo_scan.sh qsub -pe smp 1 -N HLAscan_${currentSample} -o $ANALYSISDIR/$currentSample/logs/ -e $ANALYSISDIR/$currentSample/logs/ -v BAM=$ANALYSISDIR/$currentSample/$currentSample.dedup.bam,SORTIE=$ANALYSISDIR/$currentSample/HLAscan,GeneList=$GL,LOGFILE=$ANALYSISDIR/$currentSample/logs/process_HLAscan.$logbasename.log,CONFIGFILE=$CONFIGFILE $PIPELINEBASE/haplo_scan.sh
done done
echo "End : $(date +"%F_%H-%M-%S")" echo "End : $(date +"%F_%H-%M-%S")"
varcallSNP:
# Variant calling first step # Variant calling first step
echo "### Variant calling SNP first step ###" echo "### Variant calling SNP first step ###"
echo "Start : $(date +"%F_%H-%M-%S")" echo "Start : $(date +"%F_%H-%M-%S")"
...@@ -95,6 +96,7 @@ do ...@@ -95,6 +96,7 @@ do
done done
echo "End : $(date +"%F_%H-%M-%S")" echo "End : $(date +"%F_%H-%M-%S")"
varcallCNV:
# Variant calling second step # Variant calling second step
echo "### Variant calling CNV second step ###" echo "### Variant calling CNV second step ###"
echo "Start : $(date +"%F_%H-%M-%S")" echo "Start : $(date +"%F_%H-%M-%S")"
...@@ -106,22 +108,24 @@ do ...@@ -106,22 +108,24 @@ do
done done
echo "End : $(date +"%F_%H-%M-%S")" echo "End : $(date +"%F_%H-%M-%S")"
createHTML:
# Generate a HTML report_vcf # Generate a HTML report_vcf
echo "### Generate HTML report ###" echo "### Generate HTML report ###"
echo "Start : $(date +"%F_%H-%M-%S")" echo "Start : $(date +"%F_%H-%M-%S")"
for currentSample in $samples for currentSample in $samples
do do
echo "Command : qsub -pe smp 1 -hold_jid varcallCNV_${currentSample},varcallSNP_${currentSample},HLAscan_${currentSample} -N createHTML_${currentSample} -o $ANALYSISDIR/$currentSample/logs/ -e $ANALYSISDIR/$currentSample/logs/ -v VCF_SNP=$ANALYSISDIR/$currentSample/$currentSample.pharmacoAnnot.SNP.vcf,VCF_CNV=$ANALYSISDIR/$currentSample/$currentSample.CNV.vcf,REPORTFILE=$ANALYSISDIR/$currentSample/HLAscan/Report,BAMFILE=$ANALYSISDIR/$currentSample/$currentSample.dedup.bam,TSVFILE=$TSVFILE,OUTPUTFILE=$ANALYSISDIR/$currentSample/$currentSample.html,LOGFILE=$ANALYSISDIR/$currentSample/logs/create_html.$logbasename.log,CONFIGFILE=$CONFIGFILE $PIPELINEBASE/wrapper_create_html.sh" echo "Command : qsub -pe smp 1 -hold_jid varcallCNV_${currentSample},varcallSNP_${currentSample},HLAscan_${currentSample} -N createHTML_${currentSample} -o $ANALYSISDIR/$currentSample/logs/ -e $ANALYSISDIR/$currentSample/logs/ -v VCF_SNP=$ANALYSISDIR/$currentSample/$currentSample.pharmacoAnnot.SNP.vcf,VCF_CNV=$ANALYSISDIR/$currentSample/$currentSample.CNV.vcf,REPORTFILE=$ANALYSISDIR/$currentSample/HLAscan/Report,BAMFILE=$ANALYSISDIR/$currentSample/$currentSample.dedup.bam,TSVFILE=$TSVFILE,OUTPUTFILE=$ANALYSISDIR/$currentSample/$currentSample.html,LOGFILE=$ANALYSISDIR/$currentSample/logs/create_html.$logbasename.log,CONFIGFILE=$CONFIGFILE $PIPELINEBASE/wrapper_create_html.sh"
qsub -pe smp 1 -hold_jid varcallCNV_${currentSample},varcallSNP_${currentSample},HLAscan_${currentSample} -N createHTML_${currentSample} -o $ANALYSISDIR/$currentSample/logs/ -e $ANALYSISDIR/$currentSample/logs/ -v VCF_SNP=$ANALYSISDIR/$currentSample/$currentSample.pharmacoAnnot.SNP.vcf,VCF_CNV=$ANALYSISDIR/$currentSample/$currentSample.CNV.vcf,REPORTFILE=$ANALYSISDIR/$currentSample/HLAscan/Report,BAMFILE=$ANALYSISDIR/$currentSample/$currentSample.dedup.bam,TSVFILE=$TSVFILE,OUTPUTFILE=$ANALYSISDIR/$currentSample/$currentSample.html,LOGFILE=$ANALYSISDIR/$currentSample/logs/create_html.$logbasename.log,CONFIGFILE=$CONFIGFILE $PIPELINEBASE/wrapper_create_html.sh qsub -pe smp 1 -hold_jid varcallCNV_${currentSample} -N createHTML_${currentSample} -o $ANALYSISDIR/$currentSample/logs/ -e $ANALYSISDIR/$currentSample/logs/ -v VCF_SNP=$ANALYSISDIR/$currentSample/$currentSample.pharmacoAnnot.SNP.vcf,VCF_CNV=$ANALYSISDIR/$currentSample/$currentSample.CNV.vcf,REPORTFILE=$ANALYSISDIR/$currentSample/HLAscan/Report,BAMFILE=$ANALYSISDIR/$currentSample/$currentSample.dedup.bam,TSVFILE=$TSVFILE,OUTPUTFILE=$ANALYSISDIR/$currentSample/$currentSample.html,LOGFILE=$ANALYSISDIR/$currentSample/logs/create_html.$logbasename.log,CONFIGFILE=$CONFIGFILE $PIPELINEBASE/wrapper_create_html.sh
done done
echo "End : $(date +"%F_%H-%M-%S")" echo "End : $(date +"%F_%H-%M-%S")"
MergeVCF:
#merge vcf fies #merge vcf fies
echo "### Merge VCF files" echo "### Merge VCF files"
echo "Start : $(date +"%F_%H-%M-%S")" echo "Start : $(date +"%F_%H-%M-%S")"
for currentSample in $samples for currentSample in $samples
do do
echo "Command : qsub -pe smp 1 -hold_jid varcallCNV_${currentSample},varcallSNP_${currentSample} -N mergeVCF_${currentSample} -o $ANALYSISDIR/$currentSample/logs/ -e $ANALYSISDIR/$currentSample/logs/ -v INPUTFILE1=$ANALYSISDIR/$currentSample/$currentSample.CNV.vcf,INPUTFILE2=$ANALYSISDIR/$currentSample/$currentSample.pharmacoAnnot.SNP.vcf,REF=$REF,OUTPUTFILE=$ANALYSISDIR/$currentSample/$currentSample.pharmacoAnnot.merge.vcf,LOGFILE=$ANALYSISDIR/$currentSample/logs/merge_file.$logbasename.log,CONFIGFILE=$CONFIGFILE $PIPELINEBASE/merge_GATK.sh" echo "Command : qsub -pe smp 1 -hold_jid varcallCNV_${currentSample},varcallSNP_${currentSample} -N mergeVCF_${currentSample} -o $ANALYSISDIR/$currentSample/logs/ -e $ANALYSISDIR/$currentSample/logs/ -v INPUTFILE1=$ANALYSISDIR/$currentSample/$currentSample.CNV.vcf,INPUTFILE2=$ANALYSISDIR/$currentSample/$currentSample.pharmacoAnnot.SNP.vcf,REF=$REF,OUTPUTFILE=$ANALYSISDIR/$currentSample/$currentSample.pharmacoAnnot.merge.vcf,LOGFILE=$ANALYSISDIR/$currentSample/logs/merge_file.$logbasename.log,CONFIGFILE=$CONFIGFILE $PIPELINEBASE/merge_GATK.sh"
qsub -pe smp 1 -hold_jid varcallCNV_${currentSample},varcallSNP_${currentSample} -N mergeVCF_${currentSample} -o $ANALYSISDIR/$currentSample/logs/ -e $ANALYSISDIR/$currentSample/logs/ -v INPUTFILE1=$ANALYSISDIR/$currentSample/$currentSample.CNV.vcf,INPUTFILE2=$ANALYSISDIR/$currentSample/$currentSample.pharmacoAnnot.SNP.vcf,REF=$REF,OUTPUTFILE=$ANALYSISDIR/$currentSample/$currentSample.merge.vcf,LOGFILE=$ANALYSISDIR/$currentSample/logs/merge_file.$logbasename.log,CONFIGFILE=$CONFIGFILE $PIPELINEBASE/merge_GATK.sh qsub -pe smp 1 -N mergeVCF_${currentSample} -o $ANALYSISDIR/$currentSample/logs/ -e $ANALYSISDIR/$currentSample/logs/ -v INPUTFILE1=$ANALYSISDIR/$currentSample/$currentSample.CNV.vcf,INPUTFILE2=$ANALYSISDIR/$currentSample/$currentSample.pharmacoAnnot.SNP.vcf,REF=$REF,OUTPUTFILE=$ANALYSISDIR/$currentSample/$currentSample.pharmacoAnnot.merge.vcf,LOGFILE=$ANALYSISDIR/$currentSample/logs/merge_file.$logbasename.log,CONFIGFILE=$CONFIGFILE $PIPELINEBASE/merge_GATK.sh
done done
echo "End : $(date +"%F_%H-%M-%S")" echo "End : $(date +"%F_%H-%M-%S")"
VCF_tools.py
View file @ 0208c5db
...@@ -74,5 +74,6 @@ def create_header(): ...@@ -74,5 +74,6 @@ def create_header():
##contig=<ID=chr22,length=51304566> ##contig=<ID=chr22,length=51304566>
##contig=<ID=chrX,length=155270560> ##contig=<ID=chrX,length=155270560>
##contig=<ID=chrY,length=59373566> ##contig=<ID=chrY,length=59373566>
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT dijex
""" """
return header return header
VCF_tools.pyc
View file @ 0208c5db
No preview for this file type
create_html_v2.py deleted 100755 → 0
View file @ 99b711dc
#!/usr/bin/python
### PharmAnnot PIPELINE ###
## create_html_v2.py
## Version : 1.0.0
## Description : Create a html file with vcf file, BAM file and report HLAscan
## Usage : python create_html_v2.py -v vcf_SNP -c vcf_CNV -h report_HLA -b BAM -t tsvfile -o html -e log
## Output : Annotated vcf file
## Requirements : python 2.7+
## Author : Simon.Verdez@u-bourgogne.fr
## Creation Date : 20180731
## last revision date : 20180731
## Known bugs : None
#Example : python create_html_v2.py dijex2000.vcf report_HLA BAM
import sys
import subprocess
import csv
from pysam import VariantFile
import htmltag
# Options
try:
opts, args = getopt.getopt(sys.argv[1:], 'v:h:b:o:e:')
for opt, arg in opts:
if opt in ("-v"):
vcfSNP = arg
if opt in ("-c"):
vcfCNV = arg
if opt in ("-h"):
report = arg
if opt in ("-b"):
bam = arg
if opt in ("-t"):
tsvfile = arg
elif opt in ("-o"):
sys.stdout = open(arg, 'w')
elif opt in ("-e"):
logfile = arg
except getopt.GetoptError:
print 'usage : '
sys.exit(1)
# Logging
logging.basicConfig(filename = '%s' % (logfile), filemode = 'a', level = logging.INFO, format = '%(asctime)s %(levelname)s - %(message)s')
logging.info('start')
# Manage options
if not vars().has_key('vcf_SNP'):
logging.exception('No vcf SNP file given')
sys.exit(1)
if not vars().has_key('vcf_CNV'):
logging.exception('No vcf CNV file given')
sys.exit(1)
if not vars().has_key('report'):
logging.exception('No report HLAscan given')
sys.exit(1)
if not vars().has_key('bam'):
logging.exception('No bam file given')
sys.exit(1)
if not vars().has_key('tsvfile'):
logging.exception('No tsv file given')
sys.exit(1)
# Check files
try:
f = open(vcf_SNP, 'r')
f.close()
except IOError:
logging.exception('vcf_SNP file does not exist')
sys.exit(1)
try:
f = open(vcf_CNV, 'r')
f.close()
except IOError:
logging.exception('vcf_CNV file does not exist')
sys.exit(1)
try:
f = open(report, 'r')
f.close()
except IOError:
logging.exception('Report HLAscan file does not exist')
sys.exit(1)
try:
f = open(bam, 'r')
f.close()
except IOError:
logging.exception('BAM file does not exist')
sys.exit(1)
try:
f = open(tsvfile, 'r')
f.close()
except IOError:
logging.exception('Tsv file does not exist')
sys.exit(1)
" construction du formulaire en premier "
head = htmltag.HEAD()
head = htmltag.append(htmltag.TITLE('pharmacogenetic'))
head = head.append(htmltag.STYLE(type="text/css").append("body { padding-left: 15em; font-family: Georgia; color:black }"))
head = head.append(htmltag.STYLE(type="text/css").append("h2 { padding-left: 15em;font-family: Georgia;color:blue }"))
head = head.append(htmltag.STYLE(type="text/css").append("ul.navbar { list-style-type: none; padding: 0; margin: 0; position: absolute; top: 2em; left: 1em; width: 13em }"))
head = head.append(htmltag.STYLE(type="text/css").append("ul.navbar li { background: yellow; margin: 0.5em 0; padding: 1em; border-left: 0.5em solid black }"))
head = head.append(htmltag.STYLE(type="text/css").append("P1 {color: red }"))
head = head.append(htmltag.STYLE(type="text/css").append("P2 { color: orange }"))
body = htmltag.BODY()
body = body.append(htmltag.UL(CLASS = "navbar").append(htmltag.UL(htmltag.LI(htmltag.A(href = "https://www.pharmgkb.org").append("pharmgkb")))).append(htmltag.UL(htmltag.LI(htmltag.A(href = "https://cpicpgx.org").append("cpicpgx")))).append(htmltag.UL(htmltag.LI(htmltag.A(href = "http://www.pharmacogenetics.fr/").append("pharmacogenetics.fr")))))
body = body.append(htmltag.H2('Variants pharmgkb'))
"""extraction des donnees """
data = get_data(tsvfile)
vcf_reader_SNP = VariantFile(vcfSNP)
vcf_reader_CNV = VariantFile(vcfCNV)
BAM_FILE = VariantFile(bam)
Report_HLA = open(report,"r")
for record in vcf_reader_SNP.fetch():
try:
for i in range(len(record.info["Pharmaco_Annot"])/12):
try:
gene = repr(record.info["GENEINFO"])
except:
gene = "intron"
alt = record.info["Pharmaco_Annot"][(i*12)+1].replace(")","").replace("'","").replace("-","")
pheno = record.info["Pharmaco_Annot"][(i*12)+3].replace(")","").replace("'","").replace("-"," ")
publi = record.info["Pharmaco_Annot"][(i*12)+5].replace(")","").replace("'","").replace("-"," ")
level = record.info["Pharmaco_Annot"][(i*12)+7].replace(")","").replace("'","").replace("-","")
patho = record.info["Pharmaco_Annot"][(i*12)+9].replace(")","").replace("'","").replace("-","")
molecule = record.info["Pharmaco_Annot"][(i*12)+11].replace(")","").replace("'","").replace("-"," ").replace("]","")
if level.startswith("1"):
body= body.append(htmltag.P(htmltag.P1("pathologie: "+patho))+htmltag.P(htmltag.P1("molecule: "+molecule))+htmltag.P(htmltag.P1("gene: "+gene))+htmltag.P(htmltag.P1("variant: "+repr(record.pos)+" "+repr(record.ref)+">"+repr(record.alts)+" "+repr(record.id)))+htmltag.P(htmltag.P1("base alternative: "+alt))+htmltag.P(htmltag.P1("niveau: "+level))+htmltag.P(htmltag.P1("phenotype: "+pheno)))
body= body.append(htmltag.P(htmltag.P1(htmltag.A(href = publi).append("publication: "+publi))))
elif level.startswith("2"):
body= body.append(htmltag.P(htmltag.P2("pathologie: "+patho))+htmltag.P(htmltag.P2("molecule: "+molecule))+htmltag.P(htmltag.P2("gene: "+gene))+htmltag.P(htmltag.P2("variant: "+repr(record.pos)+" "+repr(record.ref)+">"+repr(record.alts)+" "+repr(record.id)))+htmltag.P(htmltag.P2("base alternative: "+alt))+htmltag.P(htmltag.P2("niveau: "+level))+htmltag.P(htmltag.P2("phenotype: "+pheno)))
body= body.append(htmltag.P(htmltag.P2(htmltag.A(href = publi).append("publication: "+publi))))
else :
body= body.append(htmltag.P("pathologie: "+patho)+htmltag.P("molecule: "+molecule)+htmltag.P("gene: "+gene)+htmltag.P("variant: "+repr(record.pos)+" "+repr(record.ref)+">"+repr(record.alts)+" "+repr(record.id))+htmltag.P("base alternative: "+alt)+htmltag.P("niveau: "+level)+htmltag.P(htmltag.P("phenotype: "+pheno)))
body= body.append(htmltag.P(htmltag.A(href = publi).append("publication: "+publi)))
body= body.append(htmltag.P("#######################################################################"))
except:
pass
body = body.append(htmltag.H2('Deletion ou duplication'))
for record in vcf_reader_CNV.fetch():
for elm in record.alts:
if elm == "<DEL>":
ListeM = []
ListeD = []
body= body.append(htmltag.P("Deletion detectee par xHMM sur le gene "+record.info["GENE"]))
with open(data) as tsvfile:
reader = csv.DictReader(tsvfile, dialect='excel-tab')
for row in reader:
if row['gene'] != "":
if row['drug'] == record.info["GENE"]:
if row['drug'] not in ListeM:
ListeG.append(row['drug'])
for Mol in ListeM:
body= body.append(htmltag.P("######Gene ayant une relation avec :"+Mol))
depth_ave = 0.0
pos_len = 0
depth_new = 0
p = subprocess.Popen(["samtools", 'depth', BAM_FILE, '-r', record.chrom + ':' + str(record.pos) + '-'\
+ str(record.stop)], stdout = subprocess.PIPE)
for line in p.stdout.readlines():
item = line.split()
if len(item) == 3:
ListeD.append(item[2])
for elm in ListeD:
depth_new = depth_new + int(elm)
depth_ave = depth_ave + float(item[2])
pos_len = pos_len + (int(record.stop) - int(record.pos) + 1)
depth_ave = depth_new/float(pos_len)
if depth_ave > 40:
body= body.append(htmltag.P("######Detecte heterozygote, profondeur moyenne :"+depth_ave))
if depth_ave < 10:
body= body.append(htmltag.P("######Detecte homozygote, profondeur moyenne :"+depth_ave))
else:
body= body.append(htmltag.P("######Detecte probablement heterozygote, profondeur moyenne :"+depth_ave))
if elm == "<DUP>":
ListeM = []
body= body.append(htmltag.P("Duplication detectee par xHMM sur le gene "+record.info["GENE"]))
with open(data) as tsvfile:
reader = csv.DictReader(tsvfile, dialect='excel-tab')
for row in reader:
if row['gene'] != "":
if row['drug'] == record.info["GENE"]:
if row['drug'] not in ListeM:
ListeG.append(row['drug'])
for Mol in ListeM:
body= body.append(htmltag.P("######Gene ayant une relation avec :"+Mol))
HLA = Report_HLA.readlines()
ListeA = []
ListeB = []
ListeC = []
ListeDQB1 = []
ListeDRB1 = []
for lines in HLA :
if lines.startswith("A"):
for ligne in lines.split("\t"):
if ligne.startswith("A*"):
ListeA.append(ligne)
if lines.startswith("B"):
for ligne in lines.split("\t"):
if ligne.startswith("B*"):
ListeB.append(ligne)
if lines.startswith("C"):
for ligne in lines.split("\t"):
if ligne.startswith("C*"):
ListeC.append(ligne)
if lines.startswith("DQB1"):
for ligne in lines.split("\t"):
if ligne.startswith("DQB1*"):
ListeDQB1.append(ligne)
if lines.startswith("DRB1"):
for ligne in lines.split("\t"):
if ligne.startswith("DRB1*"):
ListeDRB1.append(ligne)
body = body.append(htmltag.H2('Resultat genotypage HLA'))
body = body.append(htmltag.P("HLA-A : "+ListeA[0]+" / "+ListeA[1]))
body = body.append(htmltag.P("HLA-B : "+ListeB[0]+" / "+ListeB[1]))
body = body.append(htmltag.P("HLA-C : "+ListeC[0]+" / "+ListeC[1]))
body = body.append(htmltag.P("HLA-DQB1 : "+ListeDQB1[0]+" / "+ListeDQB1[1]))
body = body.append(htmltag.P("HLA-DRB1 : "+ListeDRB1[0]+" / "+ListeDRB1[1]))
with open(data) as tsvfile:
reader = csv.DictReader(tsvfile, dialect='excel-tab')
for row in reader:
if row['gene'] != "":
if row['variants'].startswith("HLA-A"):
Adigit = row['variants'].split("*")[1]
Adigit2 = Adigit.split(":")[0]+Adigit.split(":")[1]
HL1 = ListeA[0].split("*")[1].split(":")[0]+ListeA[0].split("*")[1].split(":")[1]
HL2 = ListeA[1].split("*")[1].split(":")[0]+ListeA[1].split("*")[1].split(":")[1]
if HL1 == Adigit2:
pheno = row['effect']
publi = row['article']
level = str(row['LEVEL'])
if level.startswith("1"):
body = body.append(htmltag.P1(htmltag.P(pheno)+htmltag.P(level)))
body = body.append(htmltag.P1(htmltag.P(htmltag.A(href = publi).append("publication: "+publi))))
if level.startswith("2"):
body = body.append(htmltag.P2(htmltag.P(pheno)+htmltag.P(level)))
body = body.append(htmltag.P1(htmltag.P(htmltag.A(href = publi).append("publication: "+publi))))
else:
body = body.append(htmltag.P(pheno)+htmltag.P(level))
body = body.append(htmltag.P1(htmltag.P(htmltag.A(href = publi).append("publication: "+publi))))
if HL2 == Adigit2:
pheno = row['effect']
publi = row['article']
level = str(row['LEVEL'])
if level.startswith("1"):
body = body.append(htmltag.P1(htmltag.P(pheno)+htmltag.P(level)))
body = body.append(htmltag.P1(htmltag.P(htmltag.A(href = publi).append("publication: "+publi))))
if level.startswith("2"):
body = body.append(htmltag.P2(htmltag.P(pheno)+htmltag.P(level)))
body = body.append(htmltag.P1(htmltag.P(htmltag.A(href = publi).append("publication: "+publi))))
else:
body = body.append(htmltag.P(pheno)+htmltag.P(level))
body = body.append(htmltag.P1(htmltag.P(htmltag.A(href = publi).append("publication: "+publi))))
if row['varaints'].startswith("HLA-B"):
Bdigit = row['variants'].split("*")[1]
Bdigit2 = Bdigit.split(":")[0]+Bdigit.split(":")[1]
HL1 = ListeB[0].split("*")[1].split(":")[0]+ListeB[0].split("*")[1].split(":")[1]
HL2 = ListeB[1].split("*")[1].split(":")[0]+ListeB[1].split("*")[1].split(":")[1]
if HL1 == Bdigit2:
pheno = row['effect']
publi = row['article']
level = str(row['LEVEL'])
if level.startswith("1"):
body = body.append(htmltag.P1(htmltag.P(pheno)+htmltag.P(level)))
body = body.append(htmltag.P1(htmltag.P(htmltag.A(href = publi).append("publication: "+publi))))
if level.startswith("2"):
body = body.append(htmltag.P2(htmltag.P(pheno)+htmltag.P(level)))
body = body.append(htmltag.P2(htmltag.P(htmltag.A(href = publi).append("publication: "+publi))))
else:
body = body.append(htmltag.P(pheno)+htmltag.P(level))
body = body.append(htmltag.P(htmltag.P(htmltag.A(href = publi).append("publication: "+publi))))
if HL2 == Bdigit2:
pheno = row['effect']
publi = row['article']
level = str(row['LEVEL'])
if level.startswith("1"):
body = body.append(htmltag.P1(htmltag.P(pheno)+htmltag.P(level)))
body = body.append(htmltag.P1(htmltag.P(htmltag.A(href = publi).append("publication: "+publi))))
if level.startswith("2"):
body = body.append(htmltag.P2(htmltag.P(pheno)+htmltag.P(level)))
body = body.append(htmltag.P2(htmltag.P(htmltag.A(href = publi).append("publication: "+publi))))
else:
body = body.append(htmltag.P(pheno)+htmltag.P(level))
body = body.append(htmltag.P(htmltag.P(htmltag.A(href = publi).append("publication: "+publi))))
print htmltag.HTML(head+body)
haplo_scan.sh
View file @ 0208c5db
...@@ -15,7 +15,7 @@ ...@@ -15,7 +15,7 @@
#$ -N hlascan #$ -N hlascan
#$ -q batch #$ -q batch
#$ -pe smp 1 #$ -pe smp 2
#$ -V #$ -V
# Log file path option # Log file path option
...@@ -66,7 +66,7 @@ fi ...@@ -66,7 +66,7 @@ fi
SAMPLE=`basename ${BAM} | cut -d . -f1` SAMPLE=`basename ${BAM} | cut -d . -f1`
RES=${SORTIE} RES=${SORTIE}
GL=/work/gad/shared/analyse/test_HLA/HLAscan_v1.0/hla-ref-5gene/gene_list GL=/work/gad/shared/bin/HLAscan_v1.0/hla-ref-5gene/gene_list
#Launch haplo_scan #Launch haplo_scan
/work/gad/shared/bin/miniconda2/bin/python2.7 /work/gad/shared/bin/HLAscan_v1.0/hla-paper/haplo_scan_v4.0-sv.py $BAM $GL $RES /work/gad/shared/bin/miniconda2/bin/python2.7 /work/gad/shared/bin/HLAscan_v1.0/hla-paper/haplo_scan_v4.0-sv.py $BAM $GL $RES
......
merge_GATK.sh
View file @ 0208c5db
...@@ -79,12 +79,13 @@ then ...@@ -79,12 +79,13 @@ then
exit 1 exit 1
fi fi
JAVA=$(grep "javacmd" $CONFIGFILE | awk '{print $2}') #JAVA=$(grep "javacmd" $CONFIGFILE | awk '{print $2}')
module load java64/1.8.0_162
PATHSCRIPT=$(grep "pipelinebase" $CONFIGFILE | awk '{print $2}') PATHSCRIPT=$(grep "pipelinebase" $CONFIGFILE | awk '{print $2}')
#GATK=$(grep "GATKbase" $CONFIGFILE | awk '{print $2}') #GATK=$(grep "GATKbase" $CONFIGFILE | awk '{print $2}')
echo "command : $JAVA -jar /work/gad/shared/bin/GATK_3.8/GenomeAnalysisTK.jar -T CombineVariants -R $REF --variant $INPUTFILE1 --variant $INPUTFILE2 -o merge.vcf -genotypeMergeOptions UNIQUIFY" echo "command : $JAVA -jar /work/gad/shared/bin/GATK_3.7/GenomeAnalysisTK.jar -T CombineVariants -R $REF --variant $INPUTFILE1 --variant $INPUTFILE2 -genotypeMergeOptions UNIQUIFY"
$JAVA -jar /work/gad/shared/bin/GATK_3.8/GenomeAnalysisTK.jar -T CombineVariants -R $REF --variant $INPUTFILE2 --variant $INPUTFILE2 -genotypeMergeOptions UNIQUIFY > $OUTPUTFILE java -jar /work/gad/shared/bin/GATK_3.7/GenomeAnalysisTK.jar -T CombineVariants -R $REF --variant $INPUTFILE1 --variant $INPUTFILE2 -o $OUTPUTFILE -genotypeMergeOptions UNIQUIFY
genotypeMerge_exitcode=$? genotypeMerge_exitcode=$?
echo "GenotypeMerge exit code : $genotypeMerge_exitcode" echo "GenotypeMerge exit code : $genotypeMerge_exitcode"
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment