全长转录组 | Iso-Seq 三代测序数据分析流程 (PacBio) （3）-- SQANTI3 v5.2

$ wget https://github.com/ConesaLab/SQANTI3/archive/refs/tags/v5.2.tar.gz
$ tar -xvf v5.2.tar.gz

$ cd SQANTI3
$ conda env create -f SQANTI3.conda_env.yml
$ source activate SQANTI3.env

#将 gtfToGenePred 加入 PATH，这里根据自己的路径进行相应修改。
$ echo PATH=/mnt/data/home/mli/Desktop/Software/SQANTI3-5.2/utilities:$PATH >> ~/.bashrc && source ~/.bashrc
#将 SQANTI3 python脚本加入 PATH，这里根据自己的路径进行相应修改。
$ echo PATH=/mnt/data/home/mli/Desktop/Software/SQANTI3-5.2:$PATH >> ~/.bashrc && source ~/.bashrc

#首先激活SQANTI3
$ source activate SQANTI3.env

#然后下载cDNA_Cupcake
(SQANTI3.env)$ git clone https://github.com/Magdoll/cDNA_Cupcake.git
(SQANTI3.env)$ cd cDNA_Cupcake

#安装cDNA_Cupcake
(SQANTI3.env)$ python setup.py build
(SQANTI3.env)$ python setup.py install

(base)-bash-4.1$ conda activate SQANTI3.env
(SQANTI3.env)-bash-4.1$

#官方教程
(SQANTI3.env)-bash-4.1$ export PYTHONPATH=$PYTHONPATH:<path_to>/cDNA_Cupcake/sequence/
(SQANTI3.env)-bash-4.1$ export PYTHONPATH=$PYTHONPATH:<path_to>/cDNA_Cupcake/

#实际运行
(SQANTI3.env)-bash-4.1$ export PYTHONPATH=$PYTHONPATH:/mnt/data/home/mli/Desktop/Software/cDNA_Cupcake-master/sequence/
(SQANTI3.env)-bash-4.1$ export PYTHONPATH=$PYTHONPATH:/mnt/data/home/mli/Desktop/Software/cDNA_Cupcake-master

#进入 SQANTI3-5.2 文件夹
#修改  UHR_chr22_short_reads.fofn   中文件路径为如下：
example/UHR_Rep1_chr22.R1.fastq.gz example/UHR_Rep1_chr22.R2.fastq.gz
example/UHR_Rep2_chr22.R1.fastq.gz example/UHR_Rep2_chr22.R2.fastq.gz

#运行
$ sqanti3_qc.py example/UHR_chr22.gtf example/gencode.v38.basic_chr22.gtf example/GRCh38.p13_chr22.fasta    \
                     --CAGE_peak data/ref_TSS_annotation/human.refTSS_v3.1.hg38.bed    \
                     --polyA_motif_list data/polyA_motifs/mouse_and_human.polyA_motif.txt    \
                     -o UHR_chr22 -d example/SQANTI3_output -fl example/UHR_abundance.tsv    \
                     --short_reads example/UHR_chr22_short_reads.fofn --cpus 4 --report both  

$ python sqanti3_qc.py example/UHRR.collapsed.gff example/GRCh38.gtf example/GRCh38.fa    \
      -o UHRR -d example/SQANTI3_UHRR_Output -fl example/UHRR.collapsed.flnc_count.txt \
      --isoAnnotLite --gff3 example/human_tappas_gencode_annotation_file.gff3 \
      --cpus 20 --report both

$ conda install -c bioconda gffread
#GFF转换GTF
$ gffread input.gff3 -T -o out.gtf
#GTF转换GFF3
$ gffread input.gtf -o out.gff3

usage: sqanti3_qc.py [-h] [--min_ref_len MIN_REF_LEN] [--force_id_ignore] 
                     [--aligner_choice {minimap2,deSALT,gmap,uLTRA}] 
                     [--CAGE_peak CAGE_PEAK] [--polyA_motif_list POLYA_MOTIF_LIST]
                     [--polyA_peak POLYA_PEAK] [--phyloP_bed PHYLOP_BED] 
                     [--skipORF] [--is_fusion] [--orf_input ORF_INPUT] 
                     [--fasta] [-e EXPRESSION] [-x GMAP_INDEX] 
                     [-t CPUS] [-n CHUNKS] [-o OUTPUT] [-d DIR]
                     [-c COVERAGE] [-s SITES] [-w WINDOW] [--genename] 
                     [-fl FL_COUNT] [-v] [--saturation] 
                     [--report {html,pdf,both,skip}] 
                     [--isoAnnotLite] [--gff3 GFF3]
                     [--short_reads SHORT_READS] [--SR_bam SR_BAM] 
                     [--isoform_hits] 
                     [--ratio_TSS_metric {max,mean,median,3quartile}]
                     isoforms annotation genome

Structural and Quality Annotation of Novel Transcript Isoforms

positional arguments:
  isoforms              Isoforms (FASTA/FASTQ) or GTF format. It is recommended to 
                        provide them in GTF format, but if it is needed to map the sequences 
                        to the genome use a FASTA/FASTQ file with the
                        --fasta option.

  annotation            Reference annotation file (GTF format)

  genome                Reference genome (Fasta format)

options:
  -h, --help            show this help message and exit
  --min_ref_len MIN_REF_LEN
                        Minimum reference transcript length (default: 200 bp)

  --force_id_ignore     Allow the usage of transcript IDs non related with 
                        PacBio's nomenclature (PB.X.Y)

  --aligner_choice {minimap2,deSALT,gmap,uLTRA}

  --CAGE_peak CAGE_PEAK
                        FANTOM5 Cage Peak (BED format, optional)

  --polyA_motif_list POLYA_MOTIF_LIST
                        Ranked list of polyA motifs (text, optional)

  --polyA_peak POLYA_PEAK
                        PolyA Peak (BED format, optional)

  --phyloP_bed PHYLOP_BED
                        PhyloP BED for conservation score (BED, optional)

  --skipORF             Skip ORF prediction (to save time)

  --is_fusion           Input are fusion isoforms, must supply GTF as input

  --orf_input ORF_INPUT
                        Input fasta to run ORF on. By default, ORF is run on genome-corrected 
                        fasta - this overrides it. If input is fusion (--is_fusion), this must be provided for ORF prediction.

  --fasta               Use when running SQANTI by using as input a FASTA/FASTQ with the sequences of isoforms

  -e EXPRESSION, --expression EXPRESSION
                        Expression matrix (supported: Kallisto tsv)

  -x GMAP_INDEX, --gmap_index GMAP_INDEX
                        Path and prefix of the reference index created by gmap_build. Mandatory 
                        if using GMAP unless -g option is specified.

  -t CPUS, --cpus CPUS  Number of threads used during alignment by aligners. (default: 10)

  -n CHUNKS, --chunks CHUNKS
                        Number of chunks to split SQANTI3 analysis in for speed up (default: 1).

  -o OUTPUT, --output OUTPUT
                        Prefix for output files.

  -d DIR, --dir DIR     Directory for output files. Default: Directory where the script was run.

  -c COVERAGE, --coverage COVERAGE
                        Junction coverage files (provide a single file, comma-delmited filenames, 
                        or a file pattern, ex: "mydir/*.junctions").

  -s SITES, --sites SITES
                        Set of splice sites to be considered as canonical (comma-separated list of splice 
                        sites). Default: GTAG,GCAG,ATAC.

  -w WINDOW, --window WINDOW
                        Size of the window in the genomic DNA screened for Adenine content downstream of TTS

  --genename            Use gene_name tag from GTF to define genes. Default: gene_id used to define genes

  -fl FL_COUNT, --fl_count FL_COUNT
                        Full-length PacBio abundance file

  -v, --version         Display program version number.

  --saturation          Include saturation curves into report

  --report {html,pdf,both,skip}
                        select report format --html --pdf --both --skip

  --isoAnnotLite        Run isoAnnot Lite to output a tappAS-compatible gff3 file

  --gff3 GFF3           Precomputed tappAS species specific GFF3 file. It will serve as reference 
                        to transfer functional attributes

  --short_reads SHORT_READS
                        File Of File Names (fofn, space separated) with paths to FASTA or FASTQ from 
                        Short-Read RNA-Seq. If expression or 
                        coverage files are not provided, Kallisto (just for pair-end
                        data) and STAR, respectively, will be run to calculate them.

  --SR_bam SR_BAM       Directory or fofn file with the sorted bam files of Short Reads RNA-Seq mapped 
                        against the genome

  --isoform_hits        Report all FSM/ISM isoform hits in a separate file

  --ratio_TSS_metric {max,mean,median,3quartile}
                        Define which statistic metric should be reported in the ratio_TSS column

#教程
$ python sqanti3_filter.py rules path/to/classification.txt 

#实际运行
$ python sqanti3_filter.py rules example/SQANTI3_UHRR_Output/UHRR_classification.txt 

#教程
$ python sqanti3_filter.py ml path/to/classification.txt

#实际运行
$ python sqanti3_filter.py ml example/SQANTI3_UHRR_Output/UHRR_classification.txt 

#rules实际运行
$ python sqanti3_rescue.py rules --isoforms *_corrected.fasta --gtf *.filtered.gtf -f GRCh38.gtf -k UHRR_classification.txt

#机器学习实际运行
$ python sqanti3_rescue.py ml --isoforms *_corrected.fasta --gtf *.filtered.gtf -f GRCh38.gtf -k UHRR_classification.txt

Error compiling Cython file:
------------------------------------------------------------
...
                exon_tree.insert_interval(Interval(e_start+offset, i+offset, index))
                index += 1
                tag = False
            elif baseC[i] > 0 and (altC_pos[i] > threshSplit or altC_neg[i+1] < -threshSplit): # alt. junction found!
                # end the current exon at i and start a new one at i + 1
                print "alt. junction found at", i
                      ^
------------------------------------------------------------

cupcake/tofu/branch/c_branch.pyx:30:22: Syntax error in simple statement list
Traceback (most recent call last):
  File "/mnt/data/home/mli/Desktop/Software/cDNA_Cupcake-master/setup.py", line 25, in <module>
    ext_modules = cythonize(ext_modules),
  File "/mnt/data/home/mli/miniforge-pypy3/envs/SQANTI3.env/lib/python3.10/site-packages/Cython/Build/Dependencies.py", line 1154, in cythonize
    cythonize_one(*args)
  File "/mnt/data/home/mli/miniforge-pypy3/envs/SQANTI3.env/lib/python3.10/site-packages/Cython/Build/Dependencies.py", line 1321, in cythonize_one
    raise CompileError(None, pyx_file)
Cython.Compiler.Errors.CompileError: cupcake/tofu/branch/c_branch.pyx

编辑 setup.py文件中的 第25行：

把 ext_modules = cythonize(ext_modules),
改为 ext_modules = cythonize(ext_modules, language_level = "2"),
再次运行：
$ python setup.py build

SQANTI3.env) mli@ca496d1fda31:~$ sqanti3_qc.py 
Traceback (most recent call last):
  File "/mnt/data/home/mli/Desktop/Software/SQANTI3-5.2/sqanti3_qc.py", line 20, in <module>
    from scipy import mean
ImportError: cannot import name 'mean' from 'scipy' (/mnt/data/home/mli/miniforge-pypy3/envs/SQANTI3.env/lib/python3.10/site-packages/scipy/__init__.py)

$ conda env remove --name SQANTI3.env

$ conda env create -f SQANTI3.conda_env.yml