diff --git a/.github/workflows/build-deploy.yml b/.github/workflows/build-deploy.yml new file mode 100644 index 0000000..977b8b5 --- /dev/null +++ b/.github/workflows/build-deploy.yml @@ -0,0 +1,58 @@ +name: Create and publish a Docker image + +on: + push: + branches: main + paths-ignore: + - 'LICENSE*' + - 'README.md' + - 'requirements.txt' + - '*.yml' + release: + types: [published] + +env: + REGISTRY: ghcr.io + IMAGE_NAME: ${{ github.repository }} + +jobs: + build-and-push-image: + runs-on: ubuntu-latest + permissions: + contents: read + packages: write + + steps: + - name: Checkout repository + uses: actions/checkout@v4 + + - name: Set up Docker Buildx + uses: docker/setup-buildx-action@v3 + + - name: Log in to the Container registry + uses: docker/login-action@v3 + with: + registry: ${{ env.REGISTRY }} + username: ${{ github.actor }} + password: ${{ secrets.GITHUB_TOKEN }} + + - name: Extract metadata (tags, labels) for Docker + id: meta + uses: docker/metadata-action@v5 + with: + images: ${{ env.REGISTRY }}/${{ env.IMAGE_NAME }} + tags: | + type=raw,value=latest,enable={{is_default_branch}} + type=semver,pattern={{major}}.{{minor}}.{{patch}} + type=semver,pattern={{major}}.{{minor}} + type=semver,pattern={{major}} + + - name: Build and push Python-only Docker image + uses: docker/build-push-action@v6 + with: + context: . + push: true + tags: ${{ steps.meta.outputs.tags }} + labels: ${{ steps.meta.outputs.labels }} + cache-to: type=gha,mode=max + cache-from: type=gha diff --git a/BlueSTARR-multitask.py b/BlueSTARR-multitask.py index e786c6f..7250a8e 100755 --- a/BlueSTARR-multitask.py +++ b/BlueSTARR-multitask.py @@ -15,7 +15,7 @@ from keras.layers import Dropout, Reshape, Dense, Activation, Flatten from keras.layers import BatchNormalization, InputLayer, Input, LSTM, GRU, Bidirectional, Add, Concatenate, LayerNormalization, MultiHeadAttention import keras_nlp -from keras_nlp.layers import SinePositionEncoding +from keras_nlp.layers import SinePositionEncoding, TransformerEncoder, RotaryEmbedding from keras import models from keras.models import Sequential, Model from keras.optimizers import Adam @@ -51,7 +51,7 @@ #========================================================================= # main() #========================================================================= -def main(configFile,subdir,modelFilestem): +def main(configFile,subdir,modelFilestem, pretrained=None): startTime=time.time() #random.seed(RANDOM_SEED) @@ -62,17 +62,20 @@ def main(configFile,subdir,modelFilestem): # Load data print("loading data",flush=True) shouldRevComp=config.RevComp==1 - (X_train_sequence, X_train_seq_matrix, X_train, Y_train) = \ + (X_train_sequence, X_train_seq_matrix, X_train, Y_train, idx_train) = \ prepare_input("train",subdir,shouldRevComp,config.MaxTrain,config) - (X_valid_sequence, X_valid_seq_matrix, X_valid, Y_valid) = \ + (X_valid_sequence, X_valid_seq_matrix, X_valid, Y_valid, idx_val) = \ prepare_input("validation",subdir,shouldRevComp,config.MaxTrain,config) - (X_test_sequence, X_test_seq_matrix, X_test, Y_test) = \ + (X_test_sequence, X_test_seq_matrix, X_test, Y_test, idx_test) = \ prepare_input("test",subdir,shouldRevComp,config.MaxTest,config) \ if(config.ShouldTest!=0) else (None, None, None, None) seqlen=X_train.shape[1] # Build model model=BuildModel(seqlen) + if pretrained: + print("Loading pretrained model weights from",pretrained,flush=True) + model.load_weights(pretrained) model.summary() # Train @@ -96,7 +99,8 @@ def main(configFile,subdir,modelFilestem): numTasks=len(config.Tasks) for i in range(numTasks): summary_statistics(X_test,Y_test,"Test",i,numTasks, - config.Tasks[i],model) + config.Tasks[i],model,idx_test,modelFilestem) + print('Min validation loss:', round(min(history.history['val_loss']), 4)) # Report elapsed time endTime=time.time() @@ -106,10 +110,21 @@ def main(configFile,subdir,modelFilestem): -def summary_statistics(X, Y, set, taskNum, numTasks, taskName, model): +def summary_statistics(X, Y, set, taskNum, numTasks, taskName, model, idx, modelFilestem): pred = model.predict(X, batch_size=config.BatchSize) - cor=naiveCorrelation(Y,pred,taskNum,numTasks) + if (config.useCustomLoss) : + naiveTheta, cor=naiveCorrelation(Y,pred,taskNum,numTasks) # naiveTheta: normal scale, pred: log scale + df = pd.DataFrame({'idx':idx, 'true':tf.math.log(naiveTheta),'predicted': pred.squeeze()}) # log scale + mse = np.mean((df['true'] - df['predicted'])**2) + df.to_csv(modelFilestem+'.txt', index = False, sep='\t') + else: + cor=stats.spearmanr(tf.math.exp(pred.squeeze()),tf.math.exp(Y)) + df = pd.DataFrame({'idx':idx, 'true':Y.numpy().ravel(),'predicted': pred.squeeze()}) # log scale + mse = np.mean((df['true'] - df['predicted'])**2) + df.to_csv(modelFilestem+'.txt', index = False, sep='\t') + print(taskName+" rho=",cor.statistic,"p=",cor.pvalue) + print(taskName+' mse=', mse) @@ -118,11 +133,11 @@ def naiveCorrelation(y_true, y_pred, taskNum, numTasks): for i in range(taskNum): a+=NUM_DNA[i]+NUM_RNA[i] b=a+NUM_DNA[taskNum] c=b+NUM_RNA[taskNum] - DNA=y_true[:,a:b]+1 - RNA=y_true[:,b:c]+1 - sumX=tf.reduce_sum(DNA,axis=1) - sumY=tf.reduce_sum(RNA,axis=1) - naiveTheta=sumY/sumX + DNA=y_true[:,a:b] #+1 + RNA=y_true[:,b:c] #+1 + avgX = tf.reduce_mean(DNA, axis=1) + avgY = tf.reduce_mean(RNA, axis=1) + naiveTheta = avgY / avgX #print("naiveTheta=",naiveTheta) #print("y_pred=",tf.math.exp(y_pred[taskNum].squeeze())) cor=None @@ -130,7 +145,7 @@ def naiveCorrelation(y_true, y_pred, taskNum, numTasks): cor=stats.spearmanr(tf.math.exp(y_pred.squeeze()),naiveTheta) else: cor=stats.spearmanr(tf.math.exp(y_pred[taskNum].squeeze()),naiveTheta) - return cor + return naiveTheta, cor #def Spearman(y_true, y_pred): # return ( tf.py_function(spearmanr, [tf.cast(y_pred, tf.float32), @@ -148,10 +163,15 @@ def log(x): def logGam(x): return tf.math.lgamma(x) -def logLik(sumX,numX,Yj,logTheta,alpha,beta,numRNA): +def logLik(sumX,numX,Yj,logTheta,alpha,beta,numRNA,sumDnaLibs=None,RnaLibs=None): n=tf.shape(sumX)[0] sumX=tf.tile(tf.reshape(sumX,[n,1]),[1,numRNA]) theta=tf.math.exp(logTheta) # assume the model is predicting log(theta) + if not (sumDnaLibs is None or RnaLibs is None): + n=tf.shape(sumDnaLibs)[0] + sumDnaLibs=tf.tile(tf.reshape(sumDnaLibs,[n,1]),[1,numRNA]) + libRatio=RnaLibs/(sumDnaLibs/numX) + theta=theta*libRatio LL=(sumX+alpha)*log(beta+numX)+logGam(Yj+sumX+alpha)+Yj*log(theta)\ -logGam(sumX+alpha)-logGam(Yj+1)-(Yj+sumX+alpha)*log(theta+beta+numX) return tf.reduce_sum(LL,axis=1) # sum logLik across iid replicates @@ -159,16 +179,24 @@ def logLik(sumX,numX,Yj,logTheta,alpha,beta,numRNA): @tf.autograph.experimental.do_not_convert def makeClosure(taskNum): a=0 - for i in range(taskNum): a+=NUM_DNA[i]+NUM_RNA[i] - b=a+NUM_DNA[taskNum] - c=b+NUM_RNA[taskNum] + for i in range(taskNum): a+=NUM_DNA[i]+NUM_RNA[i] # skip previous tasks + b=a+NUM_DNA[taskNum] # first column of RNA counts + c=b+NUM_RNA[taskNum] # first column of DNA lib sizes + d=c+NUM_DNA[taskNum] # first column of RNA lib sizes @tf.autograph.experimental.do_not_convert def loss(y_true, y_pred): global EPSILON DNA=y_true[:,a:b] RNA=y_true[:,b:c] sumX=tf.reduce_sum(DNA,axis=1) - LL=-logLik(sumX,b-a,RNA,y_pred,EPSILON,EPSILON,NUM_RNA[taskNum]) + if (tf.shape(y_true)[1] > c): + DnaLibs=y_true[:,c:d] + RnaLibs=y_true[:,d:] + sumDnaLibs=tf.reduce_sum(DnaLibs,axis=1) + LL=-logLik(sumX,b-a,RNA,y_pred,EPSILON,EPSILON,NUM_RNA[taskNum], + sumDnaLibs,RnaLibs) + else: + LL=-logLik(sumX,b-a,RNA,y_pred,EPSILON,EPSILON,NUM_RNA[taskNum]) #return tf.reduce_mean(LL,axis=-1) # Put on same scale as MSE #print("pred=",y_pred) #print("-LL=",LL) @@ -191,11 +219,6 @@ def loss(y_true, y_pred): naiveTheta=sumY/sumX mse=tf.math.reduce_mean(tf.math.square(y_pred-tf.math.log(naiveTheta)), axis=1) - #mse=tf.math.reduce_mean(tf.math.square(tf.math.exp(y_pred)-naiveTheta)) - ### axis=-1) - print("log(naiveTheta)=",tf.math.log(naiveTheta)) - print("pred=",y_pred) - print("mse=",mse) return mse #cor=stats.spearmanr(tf.math.exp(y_pred[taskNum].squeeze()),naiveTheta) #return cor @@ -251,6 +274,7 @@ def loadFasta(fasta_path, as_dict=False,uppercase=False, stop_at=None, rc=generate_complementary_sequence(rec[1]) rec[1]=rec[1]+"NNNNNNNNNNNNNNNNNNNN"+rc return pd.DataFrame({'location': [e[0] for e in fastas], + 'idx': [e[0].split(' ')[0] for e in fastas], 'sequence': [e[1] for e in fastas]}) @@ -268,9 +292,11 @@ def loadCounts(filename,maxCases,config): for line in IN: if type(line) is bytes: line = line.decode("utf-8") fields=line.rstrip().split() - fields=[int(x) for x in fields] - if(config.useCustomLoss): lines.append(fields) - else: lines.append(computeNaiveTheta(fields,DNAreps,RNAreps)) + if(config.useCustomLoss): + lines.append([int(x) for x in fields]) + else: + fields=[float(x) for x in fields] # possibly normalized data + lines.append(computeNaiveTheta(fields,DNAreps,RNAreps)) linesRead+=1 if(linesRead>=maxCases): break lines=np.array(lines) @@ -283,12 +309,12 @@ def computeNaiveTheta(line,DNAreps,RNAreps): for i in range(numTasks): b=a+DNAreps[i] c=b+RNAreps[i] - DNA=line[a:b] #+1 - RNA=line[b:c] #+1 - sumX=sum(DNA)+1 - sumY=sum(RNA)+1 - naiveTheta=sumY/sumX - rec.append(naiveTheta) + DNA=line[a:b] + RNA=line[b:c] + avgX=sum(DNA)/DNAreps[i] + avgY=sum(RNA)/RNAreps[i] # normalized data + naiveTheta=avgY/avgX + rec.append(tf.math.log(naiveTheta)) # log-scale a=c return rec @@ -310,7 +336,7 @@ def prepare_input(set,subdir,shouldRevComp,maxCases,config): NUM_RNA=RNAreps matrix=pd.DataFrame(Y) matrix=tf.cast(matrix,tf.float32) - return (input_fasta_data_A.sequence, seq_matrix_A, X_reshaped, matrix) + return (input_fasta_data_A.sequence, seq_matrix_A, X_reshaped, matrix, input_fasta_data_A.idx) def BuildModel(seqlen): @@ -342,13 +368,15 @@ def BuildModel(seqlen): # Optional attention layers if(config.NumAttn>0): - x=x+keras_nlp.layers.SinePositionEncoding()(x) + x=x+keras_nlp.layers.RotaryEmbedding(name="rotary_embedding")(x) for i in range(config.NumAttn): skip=x - x=LayerNormalization()(x) - x=MultiHeadAttention(num_heads=config.AttnHeads[i], - key_dim=config.AttnKeyDim[i])(x,x) - x=Dropout(config.DropoutRate)(x) + x=LayerNormalization(name=f"layer_norm_{i}")(x) + x = TransformerEncoder(intermediate_dim=config.AttnKeyDim[i], + num_heads=config.AttnHeads[i], + dropout=config.DropoutRate, + name=f"transformer_encoder_{i}")(x) + x=Dropout(config.DropoutRate, name=f"dropout_{i}")(x) if(config.AttnResidualSkip!=0): x=Add()([x,skip]) @@ -363,7 +391,7 @@ def BuildModel(seqlen): x = Flatten()(x) # Commented out on 3/22/2023 # dense layers - if(config.NumDense>0): + if((config.NumDense>0) and config.PreDenseDropout): x=Dropout(config.DropoutRate)(x) for i in range(config.NumDense): x=kl.Dense(config.DenseSizes[i])(x) @@ -403,9 +431,12 @@ def train(model,X_train,Y_train,X_valid,Y_valid): #========================================================================= # Command Line Interface #========================================================================= -if(len(sys.argv)!=4): - exit(ProgramName.get()+" \n") -(configFile,subdir,modelFilestem)=sys.argv[1:] -main(configFile,subdir,modelFilestem) - - +import argparse +parser = argparse.ArgumentParser(prog=ProgramName.get(), description="Train a BlueSTARR model.") +parser.add_argument("configFile", help="Configuration file with hyperparameters.") +parser.add_argument("subdir", help="Subdirectory containing training, validation, and test data.") +parser.add_argument("modelFilestem", help="File (and path) stem for saving the model.") +parser.add_argument("--pretrained", help="Load pretrained model weights from this file.") +args = parser.parse_args() + +main(args.configFile,args.subdir,args.modelFilestem, pretrained=args.pretrained) \ No newline at end of file diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 0000000..6f79892 --- /dev/null +++ b/Dockerfile @@ -0,0 +1,15 @@ +ARG TENSORFLOW_VERSION=2.15.0 + +FROM tensorflow/tensorflow:${TENSORFLOW_VERSION}-gpu +ARG TENSORFLOW_VERSION + +RUN apt-get update && apt-get install --no-install-recommends -y git + +WORKDIR /bluestarr + +COPY non-tensorflow-reqs.txt . + +# We have to specify the tensorflow version or else keras_nlp will force an upgrade to the latest +# tensorflow version +RUN pip install --no-cache-dir -r non-tensorflow-reqs.txt tensorflow==${TENSORFLOW_VERSION} +COPY *.py . diff --git a/K562.config b/K562.config index aa980d9..b5e3c16 100644 --- a/K562.config +++ b/K562.config @@ -24,6 +24,7 @@ GlobalMaxPool = 0 GlobalAvePool = 1 NumDense = 0 DenseSizes = 0 +PreDenseDropout = 1 NumAttentionLayers = 0 AttentionHeads = 0 AttentionKeyDim = 0 diff --git a/NeuralConfig.py b/NeuralConfig.py index e8c27ef..271b600 100755 --- a/NeuralConfig.py +++ b/NeuralConfig.py @@ -72,6 +72,9 @@ def __init__(self,filename): self.NumDense=int(config.lookupOrDie("NumDense")) self.DenseSizes=[int(x) for x in config.lookupListOrDie("DenseSizes")] self.checkSize("DenseSizes",self.DenseSizes,self.NumDense) + self.PreDenseDropout=config.lookup("PreDenseDropout") + if(self.PreDenseDropout is None): self.PreDenseDropout=1 # defaults to true! + else: self.PreDenseDropout=int(self.PreDenseDropout) # Parameters pertaining to output and loss self.Tasks=config.lookupList("Tasks") @@ -120,12 +123,17 @@ def dump(self): print("DropoutRate=",self.DropoutRate,sep="") print("LearningRate=",self.LearningRate,sep="") print("ConvDropout=",self.ConvDropout,sep="") - + print("PreDenseDropout=",self.PreDenseDropout,sep="") + #========================================================================= # main() #========================================================================= -#if(len(sys.argv)!=2): -# exit(ProgramName.get()+" \n") -#(infile,)=sys.argv[1:] -#config=NeuralConfig(infile) -#config.dump() +def main(): + if(len(sys.argv)!=2): + exit(ProgramName.get()+" \n") + (infile,)=sys.argv[1:] + config=NeuralConfig(infile) + config.dump() + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/README.md b/README.md index a014da7..a528034 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,92 @@ # BlueSTARR BlueSTARR: predicting effects of regulatory variants + +## Running BlueSTARR + +Running the BlueSTARR code requires Python 3.10+ and the presence of a number of dependencies. These can be installed locally in the form of a conda environment, or alternatively use the Docker container automatically built and updated for the codebase. + +### Conda Environment Setup + +#### Creating a conda environment from scratch + +1. Create a new conda environment with Python 3.10. (Later versions may work but have not been tested.) + ```bash + # you can also use -n YourEnvName instead of -p /path/to/env + conda create -p /path/to/new/conda/env -c conda-forge python=3.11 + ``` +2. Activate the conda environment you just created: + ```bash + conda activate /path/to/new/conda/env + ``` +3. Install the python dependencies using `pip`: + ```bash + pip install -r requirements.txt + ``` + This assumes that you have the repository checked out and are issuing this command from the repo's root directory. + You can also do this without the full repository; all that's needed is the following two files: + - [`requirements.txt`](requirements.txt) + - [`non-tensorflow-reqs.txt`](non-tensorflow-reqs.txt) + +#### (Alternative) Cloning a previously exported conda environment + +You can “clone” (install every dependency and package at the exact same version, whether there’s a more recent compatible one or not) an existing conda environment by first exporting the configuration of an existing one (either specify it using -n/--name or -p/--prefix, or activate it first) into a file like so: +``` +conda env export > /path/to/environment.yml +``` +The [`environment.yml`](environment.yml) in this repository was created in this way. To recreate a conda environment from it, do the following: +``` +conda env create -p /path/to/new/env -f environment.yml +``` +This is not yet a complete environment because the code in this repository depends on some components of a collection of utility classes created by the Majoros lab in the past. To installl these, issue the following command (**after** you activated the conda environment you created in the previous step): +```bash +pip install "git+https://github.com/Duke-GCB/majoros-python-utils.git" +``` + +_**Caveat:** Although this used to result in a working conda environment at least with conda v4.x (which is fairly old at this point), the conda version provided by a more modern miniconda (which can also be user-installed on an HPC) apparently now fails at successfully creating the environment with this approach. Consider using the from-scratch method._ + +### Using the Docker Container + +A Docker container image is automatically built and updated from the [provided Dockerfile](Dockerfile) when the code in this repository changes or a new release is created. + +The container images are [available from the GitHub package registry](https://github.com/Duke-IGVF/BlueSTARR/pkgs/container/bluestarr). Available tags (latest, as well as major, minor and patch releases) and the command for pulling the container image can be found there. + +### Fixing TensorRT not found issue + +If your system supports TensorRT, it should be installed automatically as a dependency of `tensorflow[and-cuda]`. + +If TensorRT is installed (you can try `import tensorrt` and `import tensorrt_libs` in a Python shell; be sure to have the environment created above activated) and yet Tensorflow reports not finding it, this is likely because of a failure to load the TensorRT libraries. There are different potential causes for this: + +- The location of the libraries is not among the directories where dynamic libraries are loaded from. You can fix this by sourcing the [tensorrt-libloc.sh](tensorrt-libloc.sh) script in this repository **before** launching Python (but **after** activating your conda environment), like so: + ``` + # or use . as shorthand for source + source tensorrt-libloc.sh + ``` + +- The TensorRT libraries are not all provided with the full version (but, for example, only their major version), and Tensorflow looks for them by their full version. To fix this, run the script [tensorrt-fix.sh](tensorrt-fix.sh) (**after** activating your coda environment): + ``` + bash tensorrt-fix.sh + ``` + You need to do this only _once_ for a given conda environment. + +## Preparing Input Data for Training BlueSTARR + +The BlueSTARR model is trained on STARR-seq data, specifically DNA and RNA read counts. Typically these are obtained from the experiment in the form of [FASTQ](https://en.wikipedia.org/wiki/FASTQ_format). The FASTQ files are first converted to [BigWig](https://genome.ucsc.edu/goldenpath/help/bigWig.html) format, and then in several aggregation and filtering steps to the count tables. + +### Processing STARRseq FASTQ files to BigWig format + +To process FASTQ files into BigWig files, we use the [STARR-seq_pipeline](https://github.com/ReddyLab/cwl-pipelines/tree/main/v1.0/STARR-seq_pipeline) defined in [Common Workflow Language](https://www.commonwl.org) (CWL). + +To execute CWL workflows, you will need a CWL engine, for example [`cwltool`](https://www.commonwl.org/user_guide/introduction/prerequisites.html#cwl-runner). + +### Processing STARRseq BigWig files to count data and FASTA sequences + +For generating the counts data file and FASTA sequences for training the model, the [WGSTARR_data_preprocessing notebook](https://github.com/Duke-IGVF/BlueSTARR/edit/notebook_readme/README.md#:~:text=WGSTARR_data_preprocessing.ipynb) was used. For each 300 bp training sequence, count files were generated from the per-replicate bigWigs produced in the previous step. These files contained signal values normalized by library size, for every replicate of both the input and output libraries. + +### Downsampling + +In its current implementation, the BlueSTARR training code loads the full dataset into memory prior to commencing the model training loop. If this requires more memory than available by your compute resources, you can downsample the dataset. + +We employ two downsampling strategies that differ in the acceptance probability for each record of counts (= each training sequence): + +- **Unbiased downsampling:** each sequence is accepted with a fixed probability _N/M_, where _N_ is the desired sample size and _M_ is the total number of records being sampled. See the script [`downsample-nonuniform.py` in BlueSTARR_Evaluation_K562](https://github.com/Duke-IGVF/BlueSTARR_Evaluation_K562/blob/main/leave-one-out/BlueSTARR/leave-one-out/downsampling/downsample-nonuniform.py). +- **Biased downsampling:** the acceptance probability is a function of the estimated frequency (or kernel density) distribution of the observed activation signal $\theta$ (RNA over DNA). The script [`downsample.py` in BlueSTARR_Evaluation_K562](https://github.com/Duke-IGVF/BlueSTARR_Evaluation_K562/blob/main/leave-one-out/BlueSTARR/leave-one-out/downsampling/downsample.py) uses the empirical histogram-based PDF, where the acceptance probability for a record is $min(1, N/(M*B*p_i))$, where $p_i$ is the observed proportion of records in histogram bin _i_ into which the observed value of $\theta$ falls. The script [`downsample-biased.py` in BlueSTARR_Evaluation_A549](https://github.com/Duke-IGVF/BlueSTARR_Evaluation_A549/blob/main/full-set/BlueSTARR/downsample-biased.py) implements other functions, including PDFs and CDFs of lognormal distributions as well as powers of the lognormal CDF. Their implementation is adapted from [`Mixture-biased-sampling.ipynb` in BlueSTARR-viz](https://github.com/Duke-IGVF/BlueSTARR-viz/blob/main/sim/Mixture-biased-sampling.ipynb), where the activation-dependent acceptance probabilities and resulting enrichments in positive activations are also visualized. diff --git a/WGSTARR_data_preprocessing.ipynb b/WGSTARR_data_preprocessing.ipynb new file mode 100644 index 0000000..b2826a7 --- /dev/null +++ b/WGSTARR_data_preprocessing.ipynb @@ -0,0 +1,189 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Generate BlueSTARR training data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Whole-genome STARR-seq input (DNA, 3 reps) and output (RNA, 3 reps) data are summarized in 300bp overlapping windows (50bp step). We use the bwtool command-line tool to extract counts per window and awk to compute the average coverage per window." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%%bash\n", + "bwtool window 300 \\\n", + " ..bw \\\n", + " -step=50 \\\n", + " -skip-NA \\\n", + "| awk -vOFS=\"\\t\" 'NR==1{ split($4, vv,\",\"); nwins=length(vv)}{tot=0; split($4, vv,\",\"); for (ii=1; ii<=nwins; ii+=1){tot+=vv[ii]} print $1,$2,$3,tot/nwins}' \\\n", + "| awk '$4!=0' \\\n", + "> ..w300s50.bdg" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we find the set of windows with data for all replicates/conditions:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "bedgraphs = [\n", + "..w300s50.bdg,\n", + "..w300s50.bdg,\n", + "..w300s50.bdg\n", + "]\n", + "df = pd.read_csv(bedgraphs[0], sep='\\t', names = ['chrom','start','end','count'])\n", + "df.index = df['chrom'] + \"_\" + df['start'].astype('str') + \"_\" + df['end'].astype('str')\n", + "for ii in bedgraphs[1:]:\n", + " df_tmp = pd.read_csv(ii, sep='\\t', names = ['chrom','start','end','count'])\n", + " df_tmp.index = df['chrom'] + \"_\" + df['start'].astype('str') + \"_\" + df['end'].astype('str')\n", + " df = df.index.join(df_tmp.index, how='inner')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Use the common windows to filter the windows:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for ii in bedgraphs:\n", + " df_tmp = pd.read_csv(ii, sep='\\t', names = ['chrom','start','end','count'])\n", + " df_tmp.index = df['chrom'] + \"_\" + df['start'].astype('str') + \"_\" + df['end'].astype('str')\n", + " df_tmp.join(pd.DataFrame(index=df), how='inner').to_csv(ii.replace('.w300s50.bdg','.w300s50.in_common_win.bdg'), sep='\\t', index = False)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, merge across all samples and replicates, selecting only windows with enough average counts (100) across both input and output:\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%%bash\n", + "THRES=150;\n", + "paste -d\"\\t\" \\\n", + " .*.w300s50.in_common_win.bdg \\\n", + " .*.w300s50.in_common_win.bdg \\\n", + "|awk -vTHRES=${THRES} -vOFS=\"\\t\" 'BEGIN{print\"chrom\\tstart\\tend\\tinput_rep1\\tinput_rep2\\tinput_rep3\\toutput_rep1\\toutput_rep2\\toutput($4+$8+$12>THRES) && ($16+$20+$24>THRES){print $1,$2,$3,$4,$8,$12,$16,$20,$24}' \\\n", + "> combined.input_and_output.gt_${THRES}.bdg" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Compute a simple RNA/DNA log2fc with pseudocounts, in python:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "thres = 150\n", + "df = pd.read_csv(f'combined.input_and_output.gt_{thres}.bdg', sep='\\t')\n", + "df = df.astype(int, errors='ignore')\n", + "df.iloc[:, 3:] = df.iloc[:, 3:]/df.iloc[:, 3:].sum()*1e6\n", + "\n", + "log2FC_tmp = (np.log2((0.001+df.iloc[:, 6:].mean(axis=1)) / (0.001+df.iloc[:, 3:6].mean(axis=1))))\n", + "\n", + "df = None\n", + "df = pd.read_csv(f'combined.input_and_output.gt_{thres}.bdg', sep='\\t')\n", + "df = df.astype(int, errors='ignore')\n", + "df['log2FC'] = log2FC_tmp\n", + "df.to_csv(f'combined.input_and_output.gt_{thres}.log2FC.txt', sep='\\t', index=False, float_format='%.5f')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To add the nucleotide sequence information, we extract the windows sequences from the reference genome:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%%bash\n", + "tail -n+2 combined.input_and_output.gt_${THRES}.log2FC.txt \\\n", + "| cut -f1-3 \\\n", + "| bedtools getfasta \\\n", + " -fi GCA_000001405.15_GRCh38_full_analysis_set.fna \\\n", + " -bed stdin \\\n", + "> combined.input_and_output.gt_${THRES}.sequences.fa" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And finally, add it to the counts and log2fc table file:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%%bash\n", + "THRES=150\n", + "paste -d\"\\t\" combined.input_and_output.gt_${THRES}.log2FC.txt \\\n", + " <(cat <(echo sequence) <(cat\n", + "combined.input_and_output.gt_${THRES}.sequences.fa | awk 'NR%2==0')) \\\n", + " | gzip -c \\\n", + " > combined.input_and_output.gt_${THRES}.log2FC.sequence.txt.gz" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "BlueSTARR", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.10.13" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/environment.yml b/environment.yml new file mode 100644 index 0000000..9514c31 --- /dev/null +++ b/environment.yml @@ -0,0 +1,90 @@ +name: BlueSTARR +channels: + - defaults +dependencies: + - _libgcc_mutex=0.1 + - _openmp_mutex=5.1 + - bzip2=1.0.8 + - ca-certificates=2023.12.12 + - ld_impl_linux-64=2.38 + - libffi=3.4.4 + - libgcc-ng=11.2.0 + - libgomp=11.2.0 + - libstdcxx-ng=11.2.0 + - libuuid=1.41.5 + - ncurses=6.4 + - openssl=3.0.12 + - python=3.10.13 + - readline=8.2 + - setuptools=68.2.2 + - sqlite=3.41.2 + - tk=8.6.12 + - wheel=0.41.2 + - xz=5.4.5 + - zlib=1.2.13 + - pip: + - absl-py==2.1.0 + - astunparse==1.6.3 + - cachetools==5.3.2 + - certifi==2023.11.17 + - charset-normalizer==3.3.2 + - flatbuffers==23.5.26 + - gast==0.5.4 + - google-auth==2.27.0 + - google-auth-oauthlib==1.2.0 + - google-pasta==0.2.0 + - grpcio==1.60.0 + - h5py==3.10.0 + - idna==3.6 + - joblib==1.3.2 + - keras==2.15.0 + - keras-nlp==0.4.0 + - libclang==16.0.6 + - markdown==3.5.2 + - markupsafe==2.1.4 + - ml-dtypes==0.2.0 + - numpy==1.26.3 + - nvidia-cublas-cu12==12.2.5.6 + - nvidia-cuda-cupti-cu12==12.2.142 + - nvidia-cuda-nvcc-cu12==12.2.140 + - nvidia-cuda-nvrtc-cu12==12.2.140 + - nvidia-cuda-runtime-cu12==12.2.140 + - nvidia-cudnn-cu12==8.9.4.25 + - nvidia-cufft-cu12==11.0.8.103 + - nvidia-curand-cu12==10.3.3.141 + - nvidia-cusolver-cu12==11.5.2.141 + - nvidia-cusparse-cu12==12.1.2.141 + - nvidia-nccl-cu12==2.16.5 + - nvidia-nvjitlink-cu12==12.2.140 + - oauthlib==3.2.2 + - opt-einsum==3.3.0 + - packaging==23.2 + - pandas==2.2.0 + - pip==23.3.2 + - protobuf==4.23.4 + - pyasn1==0.5.1 + - pyasn1-modules==0.3.0 + - python-dateutil==2.8.2 + - pytz==2023.3.post1 + - requests==2.31.0 + - requests-oauthlib==1.3.1 + - rsa==4.9 + - scikit-learn==1.4.0 + - scipy==1.12.0 + - silence-tensorflow==1.2.1 + - six==1.16.0 + - support-developer==1.0.5 + - tensorboard==2.15.1 + - tensorboard-data-server==0.7.2 + - tensorflow==2.15.0.post1 + - tensorflow-estimator==2.15.0 + - tensorflow-hub==0.16.0 + - tensorflow-io-gcs-filesystem==0.35.0 + - tensorflow-text==2.15.0 + - termcolor==2.4.0 + - threadpoolctl==3.2.0 + - typing-extensions==4.9.0 + - tzdata==2023.4 + - urllib3==2.1.0 + - werkzeug==3.0.1 + - wrapt==1.14.1 diff --git a/non-tensorflow-reqs.txt b/non-tensorflow-reqs.txt new file mode 100644 index 0000000..58e57e6 --- /dev/null +++ b/non-tensorflow-reqs.txt @@ -0,0 +1,7 @@ +numpy +keras +keras-nlp==0.12.1 +pandas +scikit-learn +silence-tensorflow +git+https://github.com/Duke-GCB/majoros-python-utils.git diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..d3afe38 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,6 @@ +# There something funky with tensorflow[and-cuda] in version 2.15.0. So you +# need to use the "--extra-index-url https://pypi.nvidia.com" option +--extra-index-url https://pypi.nvidia.com +tensorflow[and-cuda]==2.15.0; sys_platform == 'linux' +tensorflow==2.15.0; sys_platform == 'darwin' # GPU isn't supported directly on macos +-r non-tensorflow-reqs.txt diff --git a/tensorrt-fix.sh b/tensorrt-fix.sh new file mode 100644 index 0000000..46debcd --- /dev/null +++ b/tensorrt-fix.sh @@ -0,0 +1,39 @@ +#!/bin/bash + +# Fix for TensorRT libraries not being found, resulting in +# Tensorflow not being able to load the TensorRT plugin. +# +# This script will create a symbolic link to the TensorRT libraries +# with the correct version number in the name. + +# Check if the TensorRT libraries are present +if python3 -c "import tensorrt_libs" &> /dev/null; then + tensorrt_lib=$(python -c "import os; import tensorrt_libs; print(os.path.dirname(tensorrt_libs.__file__))") + echo "TensorRT library found at: $tensorrt_lib" + + # Check if the full version library is present + libnvinfer=$(ls ${tensorrt_lib}/libnvinfer.so.* 2> /dev/null) + libbuilder=$(ls ${tensorrt_lib}/libnvinfer_builder_resource.so.* 2> /dev/null) + if [[ -z "$libbuilder" || -z "$libnvinfer" ]]; then + echo "Expected library or libraries not found, skipping fix." + exit 1 + fi + vinfer=${libnvinfer##*.so.} + vbuilderdiff=${libbuilder##*.so.${vinfer}} + if [[ -z $vbuilderdiff ]]; then + echo "Builder library version is the same as the main library, skipping fix." + else + for lib in $tensorrt_lib/lib*.so.${vinfer} ; do + # test whether target already exists + if [ -e $lib${vbuilderdiff} ]; then + echo "${vinfer}${vbuilderdiff} for "$(basename $lib)" exists, skipping." + else + echo "Creating symbolic link with version ${vinfer}${vbuilderdiff} for "$(basename $lib) + ln -s $lib $lib${vbuilderdiff} + fi + done + fi +else + echo "TensorRT library not found, skipping fix." + exit 0 +fi diff --git a/tensorrt-libloc.sh b/tensorrt-libloc.sh new file mode 100644 index 0000000..eb8355e --- /dev/null +++ b/tensorrt-libloc.sh @@ -0,0 +1,19 @@ +# This script must be sourced in the same shell where you will run Python +# with Tensorflow so that the environment variables are set. Like this: +# $ source tensorrt-libloc.sh + + +# Check if the TensorRT libraries are present +if python3 -c "import tensorrt_libs" &> /dev/null; then + tensorrt_lib=$(python -c "import os; import tensorrt_libs; print(os.path.dirname(tensorrt_libs.__file__))") + echo "TensorRT library found at: $tensorrt_lib" + echo "Adding TensorRT library path to LD_LIBRARY_PATH" + if [[ -z $LD_LIBRARY_PATH ]] ; then + export LD_LIBRARY_PATH=$tensorrt_lib + else + export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:$tensorrt_lib + fi +else + echo "TensorRT library not found, skipping." + exit 0 +fi