RNAseq running with umi_tools

1. install conda environment

   #conda config --set auto_activate_base false
   
   conda create --name rnaseq python=3.7
   
   #NOTE: mamba 确实快多了,以后都用 mamba❕
   #install packages
   conda activate rnaseq
   pip3 install deeptools
   pip3 install multiqc
   conda install -c bioconda stringtie subread gffread
   conda install -c conda-forge -c bioconda -c defaults -c r r-data.table r-gplots
   conda install -c conda-forge -c bioconda -c defaults -c r bioconductor-dupradar bioconductor-edger
   conda install nextflow=23.04
   
   conda install fq
   conda install -c bioconda umi_tools
   conda install -c bioconda rsem
   conda install -c bioconda salmon
   
   #conda install some tools
   #install R-packages, 
   conda install -c bioconda ucsc-bedclip
   conda install -c bioconda ucsc-bedgraphtobigwig
   conda install -c bioconda bioconductor-matrixgenerics
   #conda install -c bioconda bioconductor-deseq2
   conda install -c bioconda r-pheatmap
   conda install -c anaconda gawk
   
   conda install mamba -n base -c conda-forge
   conda config --add channels conda-forge
   mamba install -c bioconda salmon=1.10
   #salmon should be >= 1.10 since in those version salmon set `--validateMappings` as default.
   
   conda install -c bioconda trim-galore star=2.6.1d bioconductor-summarizedexperiment bioconductor-tximport bioconductor-tximeta bioconductor-deseq2
   mamba install -c bioconda samtools=1.9  
   mamba install -c conda-forge r-optparse r-vctrs=0.5.0
   conda install nextflow=23.04
   mamba install -c bioconda qualimap
   mamba install -c bioconda rseqc
   mamba install -c conda-forge openssl
   conda install -c bioconda ucsc-bedclip
   conda install -c bioconda bedtools
   conda update -c bioconda ucsc-bedclip
   #for DEBUG: bedClip: error while loading shared libraries: libssl.so.1.0.0: cannot open shared object file: No such file or directory
   conda update -c bioconda ucsc-bedgraphtobigwig
   # samtools should be >= 1.9 as only those have the option @
   #samtools sort \
   #      -@ 6 \
   #      -o HSV.d2_r1.sorted.bam \
   #      -T HSV.d2_r1.sorted \
   #      HSV.d2_r1.Aligned.out.bam
   

2. run UMItools without --umitools_dedup_stats, otherwise it cannot be finished in hamm.

   • Optimize UMItools parameters: Some parameters might influence the memory usage of UMItools. For example, you can try to reduce the number of allowed mismatches in the UMI sequence (--edit-distance-threshold). This will make the deduplication process less memory intensive but might also impact the results.

   • Use other deduplication tools: If the problem persists, you might need to use alternative tools for UMI deduplication which are less memory-intensive. Tools such as fgbio have a grouping and deduplication method similar to UMItools but have been reported to require less memory.

     #https://github.com/nf-core/rnaseq/issues/827
     #INFO for DEBUG: https://umi-tools.readthedocs.io/en/latest/faq.html
     #INFO for DEBUG: https://readthedocs.org/projects/umi-tools/downloads/pdf/stable/
     #https://github.com/CGATOxford/UMI-tools/issues/173
     # excessive dedup memory usage with output-stats #409 
     #https://github.com/CGATOxford/UMI-tools/issues/409
     #umi_tools 1.0.1
     #I am aware of previously closed issues:
     #excessive dedup memory usage #173
     #speed up stats #184
     #Running a single-end bam file with 3.13M reads and a 10bp (fully random) UMI.
     #Using --method=unique
     #There still seems to be a memory problem with --output-stats
     #Running with output-stats, memory usage climbs over 100GB and eventually crashes with "MemoryError".
     #Running without output-stats, job completes in about 3 minutes, with no problems.
     
     #TRY STANDALONE RUNNING: /usr/local/bin/python /usr/local/bin/umi_tools dedup -I HSV.d8_r1.transcriptome.sorted.bam -S HSV.d8_r1.umi_dedup.transcriptome.sorted.bam --method=unique --random-seed=100 
     #/home/jhuang/miniconda3/envs/rnaseq/bin/python /home/jhuang/miniconda3/envs/rnaseq/bin/umi_tools dedup -I star_salmon/HSV.d8_r1.sorted.bam -S HSV.d8_r1.umi_dedup.sorted.bam --output-stats HSV.d8_r1.umi_dedup.sorted --method=unique --random-seed=100
     
     #umitools dedup uses large amounts of memory and runs slowly. To speed it up it is recommended to only run it on a single chromosome, see the FAQ point number 4.
     #I suggest either making the --output-stats optional, or running a second round of deduplication on a single chromosome to generate the output stats.
     
     #--Human--
     #hamm
     /usr/local/bin/nextflow run rnaseq/main.nf --input samplesheet.csv --outdir results_GRCh38 --genome GRCh38   --with_umi --umitools_extract_method "regex" --umitools_bc_pattern "^(?P<umi_1>.{12}).*"         -profile docker     -resume  --max_cpus 54 --max_memory 120.GB --max_time 2400.h    --save_align_intermeds --save_unaligned --save_reference    --aligner 'star_salmon' --pseudo_aligner 'salmon'    --umitools_grouping_method 'unique'
     
     #sage
     nextflow run rnaseq/main.nf                --input samplesheet.csv --outdir results_GRCh38 --genome GRCh38   --with_umi --umitools_extract_method "regex" --umitools_bc_pattern "^(?P<umi_1>.{12}).*"         -profile test_full  -resume                --max_memory 256.GB --max_time 2400.h    --save_align_intermeds --save_unaligned --save_reference    --aligner 'star_salmon' --pseudo_aligner 'salmon'
     
     #--Virus--
     /usr/local/bin/nextflow run rnaseq/main.nf --input samplesheet.csv --outdir results_virus    --fasta "/home/jhuang/DATA/Data_Manja_RNAseq_Organoids_Virus/X14112.1.fasta" --gtf "/home/jhuang/DATA/Data_Manja_RNAseq_Organoids_Virus/X14112.1_v4.gtf"   --with_umi --umitools_extract_method "regex" --umitools_bc_pattern "^(?P<umi_1>.{12}).*" --umitools_dedup_stats    --skip_rseqc --skip_dupradar --skip_preseq    -profile test_full -resume  --max_cpus 55 --max_memory 120.GB --max_time 2400.h    --save_align_intermeds --save_unaligned --save_reference    --aligner 'hisat2'    --gtf_extra_attributes 'gene_name' --gtf_group_features 'gene_id' --featurecounts_group_type 'gene_name' --featurecounts_feature_type 'exon'    --umitools_grouping_method 'unique'
     

3. R-code for evaluation of nextflow outputs

   # Import the required libraries
   library("AnnotationDbi")
   library("clusterProfiler")
   library("ReactomePA")
   library(gplots)
   
   library(tximport)
   library(DESeq2)
   
   setwd("~/DATA/Data_Manja_RNAseq_Organoids/results_GRCh38_unique/star_salmon")
   
   # Define paths to your Salmon output quantification files
   files <- c("control_r1" = "./control_r1/quant.sf",
             "control_r2" = "./control_r2/quant.sf",
             "HSV.d2_r1" = "./HSV.d2_r1/quant.sf",
             "HSV.d2_r2" = "./HSV.d2_r2/quant.sf",
             "HSV.d4_r1" = "./HSV.d4_r1/quant.sf",
             "HSV.d4_r2" = "./HSV.d4_r2/quant.sf",
             "HSV.d6_r1" = "./HSV.d6_r1/quant.sf",
             "HSV.d6_r2" = "./HSV.d6_r2/quant.sf",
             "HSV.d8_r1" = "./HSV.d8_r1/quant.sf",
             "HSV.d8_r2" = "./HSV.d8_r2/quant.sf")
   
   # Import the transcript abundance data with tximport
   txi <- tximport(files, type = "salmon", txIn = TRUE, txOut = TRUE)
   
   # Define the replicates and condition of the samples
   replicate <- factor(c("r1", "r2", "r1", "r2", "r1", "r2", "r1", "r2", "r1", "r2"))
   condition <- factor(c("control", "control", "HSV.d2", "HSV.d2", "HSV.d4", "HSV.d4", "HSV.d6", "HSV.d6", "HSV.d8", "HSV.d8"))
   
   # Define the colData for DESeq2
   colData <- data.frame(condition=condition, replicate=replicate, row.names=names(files))
   
   # Create DESeqDataSet object
   dds <- DESeqDataSetFromTximport(txi, colData=colData, design=~condition)
   
   # In the context of your new code which is using tximport and DESeq2, you don't necessarily need this step. The reason is that DESeq2 performs its own filtering of low-count genes during the normalization and differential expression steps.
   # Filter data to retain only genes with more than 2 counts > 3 across all samples
   # dds <- dds[rowSums(counts(dds) > 3) > 2, ]
   
   # Run DESeq2
   dds <- DESeq(dds)
   
   # Perform rlog transformation
   rld <- rlogTransformation(dds)
   
   # Output raw count data to a CSV file
   write.csv(counts(dds), file="transcript_counts.csv")
   
   # -- gene-level count data --
   # Read in the tx2gene map from salmon_tx2gene.tsv
   #tx2gene <- read.csv("salmon_tx2gene.tsv", sep="\t", header=FALSE)
   tx2gene <- read.table("salmon_tx2gene.tsv", header=FALSE, stringsAsFactors=FALSE)
   
   # Set the column names
   colnames(tx2gene) <- c("transcript_id", "gene_id", "gene_name")
   
   # Remove the gene_name column if not needed
   tx2gene <- tx2gene[,1:2]
   
   # Import and summarize the Salmon data with tximport
   txi <- tximport(files, type = "salmon", tx2gene = tx2gene, txOut = FALSE)
   
   # Continue with the DESeq2 workflow as before...
   colData <- data.frame(condition=condition, replicate=replicate, row.names=names(files))
   dds <- DESeqDataSetFromTximport(txi, colData=colData, design=~condition)
   #dds <- dds[rowSums(counts(dds) > 3) > 2, ]    #60605-->26543
   dds <- DESeq(dds)
   rld <- rlogTransformation(dds)
   write.csv(counts(dds, normalized=FALSE), file="gene_counts.csv")
   
   #TODO: why a lot of reads were removed due to the too_short?
   #STAR --runThreadN 4 --genomeDir /path/to/GenomeDir --readFilesIn /path/to/read1.fastq /path/to/read2.fastq --outFilterMatchNmin 50 --outSAMtype BAM SortedByCoordinate --outFileNamePrefix /path/to/output
   
   dim(counts(dds))
   head(counts(dds), 10)
   

like unlike

点赞本文的读者

还没有人对此文章表态

本文有评论

没有评论

看文章，发评论，不要沉默