Nanopore可变剪切分析

Nanopore native RNA sequencing of a human poly(A) transcriptome
FLAIR原始文章

FLAIR分析流程示意图

Nanopore long reads测序相对与传统二代测序short reads 的测序而言，对于可变剪切的分析具有巨大的优势。
简单来说该流程可以应用short reads RNA-seq对Nanopore long reads 的测序的结果进行校正，然后可以获得相对较为准确的各个isoform的表达信息，并可以对于该结果进行统计检验。
我使用的是SLURM 递交系统的HPC，我的的超算上还没有安装该module所以我使用的是Conda管理这个软件。
PC和服务器conda的安装
我使用的是shell来运行，推荐大家使用snakemake workflow，后期改过会更新
1、根据软件说明设置python环境和安装包

使用conda管理文件提前创建flair_env环境（使用python3.6）,
source /your_conda_dir/conda/etc/profile.d/conda.sh
python 
#Users can run FLAIR within the conda environment provided in `misc/flair_conda_env.yaml`. FLAIR should run smoothly in this environment. Refer to the conda docs for how to [create an environment from an environment.yml file](https://docs.conda.io/projects/conda/en/latest/user-guide/tasks/manage-environments.html#creating-an-environment-from-an-environment-yml-file).
#misc/flair_conda_env.yaml 下载下来，用conda安装
conda env create -f ./flair_conda_env.yaml

2、将对应RNA-seq（short reads）bam文件转化为bed12格式文件

#!/bin/bash
#SBATCH --time=1:00:00
#SBATCH --cpus-per-task=2
#SBATCH --mem=2g


dir=/your_data_dir/Nanopores
data_dir=${dir}/MINIMAP2_bam
log_dir=${dir}/log/Bam2Bed12
bam2Bed12_dir=${dir}/Bam2Bed12
job_dir=${dir}/job/Bam2Bed12

mkdir -p ${log_dir}
mkdir -p ${job_dir}
mkdir -p ${bam2Bed12_dir}
script_dir=/flair_dir/flair-master/bin/bam2Bed12.py
thread=5
cd ${data_dir}
for i in $(ls $pwd *.bam);do
job_file="${job_dir}/${i%.bam*}.job"

    echo "#!/bin/bash
#SBATCH --job-name=${i%.bam*}.bam2Bed12.job
#SBATCH --output=$log_dir/${i%.bam*}.bam2Bed12.out
#SBATCH --time=5:00:00
#SBATCH --cpus-per-task=${thread}
#SBATCH --mem=10g
source /conda_dir/conda/etc/profile.d/conda.sh
conda activate flair_env
export PATH=\$PATH:/flair_dir/flair-master/bin/
flair_dir=/flair_dir/flair-master/flair.py
ml bedtools
ml samtools
ml minimap2

python ${script_dir} -i ${data_dir}/${i} >${bam2Bed12_dir}/${i%.bam*}.bed12

" > $job_file
sbatch $job_file
done

3、Short-read junctions 从RNA-seq中获取junctions数据

#!/bin/bash
#SBATCH --time=1:00:00
#SBATCH --cpus-per-task=2
#SBATCH --mem=2g


dir=your_data_dir
data_dir=${dir}/short_reads/human
log_dir=${dir}/log/junctions_from_sam
junctions_from_sam_dir=${dir}/junctions_from_sam/human
job_dir=${dir}/job/junctions_from_sam

mkdir -p ${log_dir}
mkdir -p ${job_dir}
mkdir -p ${junctions_from_sam_dir}
script_dir=/flair_dir/flair-master//bin/junctions_from_sam.py
thread=5
cd ${data_dir}
for i in $(ls $pwd *.bam);do
job_file="${job_dir}/${i%.bam*}.job"

    echo "#!/bin/bash
#SBATCH --job-name=${i%.bam*}.junctions_from_sam.job
#SBATCH --cpus-per-task=${thread}
#SBATCH --output=${log_dir}/${i%.bam*}.junctions_from_sam.out
#SBATCH --mem=20g
#SBATCH --time=20:00:00
source /conda_dir/conda/etc/profile.d/conda.sh
conda activate flair_env
export PATH=\$PATH:/flair_dir/flair-master//bin
ml minimap2
ml bedtools
ml samtools
cd ${junctions_from_sam_dir}
python ${script_dir} -s ${data_dir}/${i} -n ${i%.bam*}
" > $job_file
sbatch $job_file
done

4、flair nanopore reads correct

#!/bin/bash
#SBATCH --time=1:00:00
#SBATCH --cpus-per-task=2
#SBATCH --mem=2g


dir=/your_data_dir
genome=/ori_genomic/hg38/fasta/GRCh38.p13.genome.fa
flair_dir=/flair_dir/flair-master//flair-master/flair.py
data_dir=${dir}/Bam2Bed12/human
log_dir=${dir}/log/Correct
correct_dir=${dir}/Correct/human
job_dir=${dir}/job/Correct
short_reads_junctions_dir=/short_reads/junctions_from_sam/human

mkdir -p ${log_dir}
mkdir -p ${job_dir}
mkdir -p ${correct_dir}

thread=5
cd ${data_dir}
for i in $(ls $pwd *bed12);do
job_file="${job_dir}/${i%.minimap2.human.bed12*}.human.job"

    echo "#!/bin/bash
#SBATCH --job-name=${i%.minimap2.human.bed12*}.human.correct.job
#SBATCH --output=$log_dir/${i%.minimap2.human.bed12*}.human.correct.out
#SBATCH --time=5:00:00
#SBATCH --cpus-per-task=${thread}
#SBATCH --mem=10g
source /conda_dir/conda/etc/profile.d/conda.sh
conda activate flair_env
export PATH=\$PATH:/flair_dir/flair-master//flair-master/bin

ml bedtools
ml samtools
ml minimap2
cd ${correct_dir}
python ${flair_dir} correct -t ${thread} -q ${data_dir}/${i} -g ${genome} -o ${i%.minimap2.human.bed12*}.human \
    -j ${short_reads_junctions_dir}/${i%.minimap2.human.bed12*}.human_junctions.bed
" > $job_file
sbatch $job_file
done

collapse
If there are multiple samples to be compared, the flair-corrected read psl files should be concatenated prior to running flair-collapse. In addition, all raw read fastq/fasta files should either be specified after -r with space/comma separators or concatenated into a single file.

#!/bin/bash
#SBATCH --job-name=human.Collapse.job
#SBATCH --output=/your_data_dir/log/Collapse/human.Collapse.out
#SBATCH --time=24:00:00
#SBATCH --cpus-per-task=5
#SBATCH --mem=20g
source /your_conda_Dir/conda/etc/profile.d/conda.sh
conda activate flair_env
export PATH=$PATH:/flair_dir/flair-master//flair-master/bin

ml bedtools
ml samtools
ml minimap2

cd /your_data_dir/Collapse/human
python /flair_dir/flair-master//flair-master/flair.py collapse -t 5 -r /your_data_dir/Fastq/concatenated/All.fastq --temp_dir /your_data_dir/Collapse/tmp/human -q /your_data_dir/Correct/human/concatenated/human_all_corrected.bed -g /ori_genomic/hg38/fasta/GRCh38.p13.genome.fa -f  /ori_genomic/hg38/gtf/gencode.v33.annotation.gtf -o All.human

6.quantify

#!/bin/bash
#SBATCH --job-name=human.Collapse.job
#SBATCH --output=/your_data_dir/log/Collapse/human.Quantify.out
#SBATCH --time=24:00:00
#SBATCH --cpus-per-task=5
#SBATCH --mem=20g
source /your_conda_dir/conda/etc/profile.d/conda.sh
conda activate flair
export PATH=$PATH:/your_flair_dir/flair-master/bin

ml bedtools
ml samtools
ml minimap2

cd /your_data_dir/Quantify/human
python /your_flair_dir/flair-master/flair.py quantify -t 5 \
-r /your_data_dir/Quantify/human/reads_manifest.txt \
-i /your_data_dir/Collapse/human/All.human.isoforms.fa \
--temp_dir /your_data_dir/Collapse/tmp/human --tpm \
-o Nanopore.human

7.diffExp analysis

#!/bin/bash
#SBATCH --time=10:00:00
#SBATCH --cpus-per-task=5
#SBATCH --mem=10g
source /your_conda_dir/conda/etc/profile.d/conda.sh
conda activate flair_env
export PATH=$PATH:/your_flair_dir/flair-master/bin
flair_dir=/your_flair_dir/flair-master/flair.py
cd /your_data_dir/AS_diff/human/
python ${flair_dir} diffExp -t 5 -q /your_data_dir/Quantify/human/Nanopore.human -o diffExp

8.终于到了最后diffSplice分析

#!/bin/bash
#SBATCH --time=2:00:00
#SBATCH --cpus-per-task=5
#SBATCH --mem=10g
source /your_conda_dir/conda/etc/profile.d/conda.sh
conda activate flair_env

flair_dir=/your_flair_dir/flair-master/flair.py
cd /your_data_dir/AS_diff/human/diffSplice
python ${flair_dir} diffSplice -t 5 -q /your_data_dir/Quantify/human/Nanopore.human \
-i /your_data_dir/Collapse/human/All.human.isoforms.bed \
--test --conditionA Fed --conditionB Fasting \
-o Human_diffSplice

9.可视化

#! /bin/bash
#SBATCH --time=10:00:00
#SBATCH --cpus-per-task=5
#SBATCH --mem=10g
data_dir=/your_data_dir/AS_diff/human/plot_isoform
out_dir=${data_dir}/DiffAS
mkdir -p ${out_dir}
source /your_conda_dir/conda/etc/profile.d/conda.sh
conda activate flair_env
export PATH=$PATH:/your_flair_dir/flair-master/bin
plot_isoform_usage=/your_flair_dir/flair-master/bin/plot_isoform_usage.py

cd ${out_dir}
cat ${data_dir}/list/degs_ensg_id.tsv | while read line
do
python ${plot_isoform_usage} /your_data_dir/Collapse/human/All.human.isoforms.bed \
 /your_data_dir/Quantify/human/Nanopore.human \
${line}
done