Skip to content

Commit

Permalink
added scripts to use metaphlan vith VSCs and profile metagenomes
Browse files Browse the repository at this point in the history
  • Loading branch information
azufre451 committed Mar 1, 2021
1 parent 34c0de8 commit 52aea78
Show file tree
Hide file tree
Showing 4 changed files with 243 additions and 11 deletions.
28 changes: 17 additions & 11 deletions bt2_against_metagenomes/parse_bedtools_output_breadth_of_cov.py
Original file line number Diff line number Diff line change
@@ -1,12 +1,13 @@
#!/usr/bin/env python3

import sys
import gzip
import os
import pandas as pd
import numpy as np
from collections import Counter
res=[]
BREADTH_THR=0.5
BREADTH_THR=0.8


### ###
Expand All @@ -22,7 +23,7 @@
parser = argparse.ArgumentParser()

parser.add_argument("--unite")
parser.add_argument("--unitegrp", default='VC1',help='VC1 or MCL')
parser.add_argument("--unitegrp", default='MCL',help='VC1 or MCL')
parser.add_argument("--replace", action="store_true")
parser.add_argument("--quiet", action="store_true")
parser.add_argument("--limit", type=int)
Expand Down Expand Up @@ -55,23 +56,25 @@

for u in args.input:

target='pieces/'+os.path.basename(u.strip().replace('.csv','.pd').replace('.fcsv','.pd'))
target='pieces/'+os.path.basename(u.strip().replace('.fcsv.gz','.pd'))


if (not os.path.isfile(target) or os.stat(target).st_size == 0) or args.replace:

ii+=1

filename=u.strip()
print( os.path.basename(filename.replace('.csv','').replace('.fcsv','')) )
dataset,sample=os.path.basename(filename.replace('.csv','').replace('.fcsv','')).split('__')
print( os.path.basename(filename.replace('.fcsv.gz','')) )
print(os.path.basename(filename.replace('.fcsv.gz','')))
print(os.path.basename(filename.replace('.fcsv.gz','')).split('__'))
dataset,sample=os.path.basename(filename.replace('.fcsv.gz','')).split('__')

if not args.quiet:
print(dataset,sample,ii)


fline=0
a1=open(filename,'r')
a1=gzip.open(filename,'rt')

DSCT={}
itte=0
Expand Down Expand Up @@ -117,10 +120,10 @@



a[['dataset','sampleID','VC1','VC2','VC3','clusterType','complete_cluster','length','breadth','depth']].to_csv('pieces/'+os.path.basename(filename.replace('.csv','.pd').replace('.fcsv','.pd')),sep='\t')
a[['dataset','sampleID','VC1','VC2','VC3','clusterType','complete_cluster','length','breadth','depth']].to_csv('pieces/'+os.path.basename(filename.replace('.fcsv.gz','.pd')),sep='\t')
#else:
#print(target, 'already exists')
#sys.exit(0)
sys.exit(0)


else:
Expand Down Expand Up @@ -149,7 +152,9 @@
if args.unite == 'median':
aggf=np.median
if args.unite == 'avg':
aggf=np.mea
aggf=np.mean
if args.unite == 'max':
aggf=np.max

else:
aggf=np.max
Expand Down Expand Up @@ -181,6 +186,7 @@ def tak1(a):
#a.to_csv('all_breadth.csv',sep='\t')
#sys.exit(0)
a=a[a['body_site']=='stool']
a=a[a['breadth'] > BREADTH_THR]
a['VC1f'] = a['VC2'].apply(tak1)

if args.unitegrp == 'MCL':
Expand All @@ -203,14 +209,14 @@ def tak1(a):



a.to_csv('./vcta.csv',sep='\t')
#a.to_csv('./vcta.csv',sep='\t')
print(a.shape)
#a[a['body_site'].isnull()].to_csv('elab/lonley.csv',sep='\t')



vct1=pd.pivot_table(a,columns=['sampleID','dataset','body_site','country','non_westernized'],index=[grpField],values=['breadth','length'],aggfunc=aggf)['breadth']
vct1.fillna(0).to_csv('./all_vct1.csv',sep='\t')
vct1.fillna(0).to_csv('./all_vct1_breadth.csv',sep='\t')

#print(vct1)
#sys.exit(0)
Expand Down
46 changes: 46 additions & 0 deletions metaphlan-vir_metagenomes/launch_vir.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,46 @@
#!/bin/bash

#PBS -l place=free
#PBS -V

DBF=/shares/CIBIO-Storage/CM/scratch/users/moreno.zolfo/virome_data/mp3_release/VSC5_MP3/database/release_01

if [ ! $ncores ]; then ncores="2"; fi;

ldn=${dn//"/"/"__"};
NODE_TEMPFOLDER=/tmp/
SERVER_TEMPFOLDER=/shares/CIBIO-Storage/CM/tmp/mzolfo/tmp_data/
SERVER_BASEFOLDER=/shares/CIBIO-Storage/CM/scratch/users/moreno.zolfo/virome_data/prevalence4/${ldn}/

FREESPACE=$(df --output=avail $NODE_TEMPFOLDER | tail -n1);
WORKSIZE=$(du -L -ck ${prefix}/${dn}/reads/${fn}/*${extension} | tail -n1 | cut -f1)
REQSIZE=$((WORKSIZE*6))

if [[ $FREESPACE -lt $REQSIZE ]]; then
TEMPFOLDER=${SERVER_TEMPFOLDER}/MP3_${ldn}___${fn};
else
TEMPFOLDER=${NODE_TEMPFOLDER}/MP3_${ldn}___${fn};
fi

if [ -d ${TEMPFOLDER} ]; then rm -r ${TEMPFOLDER}/; fi;

mkdir -p $TEMPFOLDER/
mkdir -p $SERVER_BASEFOLDER/

$uncompress_cmd ${prefix}/${dn}/reads/${fn}/*${extension} | metaphlan --input_type fastq --bowtie2db ${DBF} -x mpa_v30_CHOCOPhlAn_201901_vsc --profile_vsc --vsc_out ${TEMPFOLDER}/${fn}.vsc.tsv --nproc ${ncores} --bowtie2out ${TEMPFOLDER}/${fn}.bowtie.bz2 > ${TEMPFOLDER}/${fn}.profile.tsv

mv ${TEMPFOLDER}/*.tsv ${SERVER_BASEFOLDER}/;
rm ${SERVER_BASEFOLDER}/${dn}__${fn}.wk
rm -r ${TEMPFOLDER};












109 changes: 109 additions & 0 deletions metaphlan-vir_metagenomes/map_launch.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,109 @@
#!/bin/bash
curDir=$(realpath $(dirname $0));

source ${curDir}/workunit_data.sh

##############################################

mode=$1;

if [ $mode == "U" ]; then
extension='.fastq'
extraction_cmd='cat'
fi;

if [ $mode == "bzip" ]; then
extension='.fastq.bz2'
extraction_cmd='bzcat'
fi;

if [ $mode == "gzip" ]; then
extension='.fastq.gz'
extraction_cmd='zcat'
fi;

if [ $mode == "fqgz" ]; then
extension='.fq.gz'
extraction_cmd='zcat'
fi;


for et in $etp; do
for utp in $(find ${prefix}/${et}/reads -maxdepth 1 -mindepth 1 -type d); do
e+=($utp);
done;

mkdir -p ${logFolder}/${et}/;
mkdir -p ${odir}/${et}/

done;

while true; do

bt=0;
b_short=$(qstat -u moreno.zolfo | grep short | grep "VDBP" | wc -l);
b_common=$(qstat -u moreno.zolfo | grep common | grep "VDBP" | wc -l);
b_cibio=$(qstat -u moreno.zolfo | grep CIBIO_ | grep "VDBP" | wc -l);
b_cibiocm=$(qstat -u moreno.zolfo | grep CIBIOCM | grep "VDBP" | wc -l);


for folder in ${e[@]}; do


sample=$(basename $folder);
dataset=$(basename $(dirname $(dirname $folder)));
odirE=${odir}/${dataset}/


if [ ! -f ${odir}/${dataset}/${sample}.vsc.tsv ] && [ ! -f ${odir}/${dataset}/${dataset}__${sample}.wk ]; then

if [ $b_short -lt 30 ]; then

echo "QUEUE (short)" ${dataset} ${sample};
qsub -q short_cpuQ -v prefix=\"${base}\",dn=\"${dataset}\",fn=\"${sample}\",uncompress_cmd=\"${extraction_cmd}\",extension=\"${extension}\",ncores=\"8\" -N VDBP_${dataset}_${sample} -o ${logFolder}/${dataset}/VDBP_${dataset}_${sample}.o.log -e ${logFolder}/${dataset}/VDBP_${dataset}_${sample}.e.log -l select=1:ncpus=8 ${curDir}/launch_vir.sh;
touch ${odir}/${dataset}/${dataset}__${sample}.wk
b_short=$((b_short+1));
bt=$((bt+1));

elif [ $b_common -lt 60 ]; then

echo "QUEUE (common)" ${dataset} ${sample};
qsub -q common_cpuQ -v prefix=\"${base}\",dn=\"${dataset}\",fn=\"${sample}\",uncompress_cmd=\"${extraction_cmd}\",extension=\"${extension}\",ncores=\"2\" -N VDBP_${dataset}_${sample} -o ${logFolder}/${dataset}/VDBP_${dataset}_${sample}.o.log -e ${logFolder}/${dataset}/VDBP_${dataset}_${sample}.e.log -l select=1:ncpus=2 ${curDir}/launch_vir.sh;
touch ${odir}/${dataset}/${dataset}__${sample}.wk
b_common=$((b_common+1));
bt=$((bt+1));

elif [ $b_cibio -lt 25 ]; then

echo "QUEUE (common)" ${dataset} ${sample};
qsub -q CIBIO_cpuQ -v prefix=\"${base}\",dn=\"${dataset}\",fn=\"${sample}\",uncompress_cmd=\"${extraction_cmd}\",extension=\"${extension}\",ncores=\"2\" -N VDBP_${dataset}_${sample} -o ${logFolder}/${dataset}/VDBP_${dataset}_${sample}.o.log -e ${logFolder}/${dataset}/VDBP_${dataset}_${sample}.e.log -l select=1:ncpus=2 ${curDir}/launch_vir.sh;
touch ${odir}/${dataset}/${dataset}__${sample}.wk
b_cibio=$((b_cibio+1));
bt=$((bt+1));


elif [ $b_cibiocm -lt 25 ]; then

echo "QUEUE (common)" ${dataset} ${sample};
qsub -q CIBIOCM_cpuQ -v prefix=\"${base}\",dn=\"${dataset}\",fn=\"${sample}\",uncompress_cmd=\"${extraction_cmd}\",extension=\"${extension}\",ncores=\"2\" -N VDBP_${dataset}_${sample} -o ${logFolder}/${dataset}/VDBP_${dataset}_${sample}.o.log -e ${logFolder}/${dataset}/VDBP_${dataset}_${sample}.e.log -l select=1:ncpus=2 ${curDir}/launch_vir.sh;
touch ${odir}/${dataset}/${dataset}__${sample}.wk
b_cibiocm=$((b_cibiocm+1));
bt=$((bt+1));

fi;
# else
# echo "SKIP" ${dataset} ${sample};
fi;

done;

a=$(qstat -u moreno.zolfo | grep "VDBP" | grep " R " | wc -l);
q=$(qstat -u moreno.zolfo | grep "VDBP" | grep " Q " | wc -l);
ete=${#e[@]}
ft=$(find ${odir}/ -name "*.vsc.tsv" | wc -l);

tg=#$(($ete-$ft));
echo $(date) "| " $bt "jobs launched, " $a "running ", $q "queued " $tg" to go. See you in one minute." ;
sleep 60
b=0;
done;
71 changes: 71 additions & 0 deletions metaphlan-vir_metagenomes/workunit_data.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,71 @@
#!/bin/bash

odir="/shares/CIBIO-Storage/CM/scratch/users/moreno.zolfo/virome_data/prevalence4/";
logFolder="/shares/CIBIO-Storage/CM/scratch/users/moreno.zolfo/virome_data/prevalence4/scripts/logs";
prefix="/shares/CIBIO-Storage/CM/cmstore/data/meta"
base=$prefix;

etp='HMP_2012
LiuB_2012
BengtssonPalmeJ_2015
PiperHG_2016
ContevilleLC_2019
GopalakrishnanV_2018
RosaBA_2018
KieserS_2018
LassalleF_2017
SankaranarayananK_2015
RampelliS_2015
MatsonV_2018
SmitsSA_2017
WingleeK_2017
LomanNJ_2013
OlmMR_2017
LawrenceA_2015
GeversD_2014
Heitz-BuschartA_2016
LiSS_2016
LokmerA_2019
Obregon-TitoAJ_2015
GuptaA_2019
WampachL_2018
VoigtAY_2015
RaymondF_2016
YeZ_2018
KorpelaK_2016
ThomasAM_2019_c
HanniganGD_2017
LoombaR_2017
LouisS_2016
ChengpingW_2017
IjazUZ_2017
ChuDM_2017
DhakanDB_2019
LiuW_2016
VogtmannE_2016
HeQ_2017
KosticAD_2015
WirbelJ_2018
YuJ_2015
KarlssonFH_2013
FengQ_2015
ZellerG_2014
WenC_2017
LiJ_2017
HansenLBS_2018
VincentC_2016
QinN_2014
XieH_2016
LiJ_2014
HallAB_2017
CosteaPI_2017
PehrssonE_2016
YassourM_2018
LeChatelierE_2013
QinJ_2012
JieZ_2017
NielsenHB_2014
BackhedF_2015
SchirmerM_2016
VatanenT_2016
YachidaS_2019';

0 comments on commit 52aea78

Please sign in to comment.