首先需要自己获取样本的SRR号,放在一个文件夹里
示例:
fanhz1@yanglab:~/experiment/Arabidopsis_thaliana$ cat Arabidopsis_thaliana.txt
SRR13711353
SRR13711354
SRR13711356
SRR13711357
SRR13711359
SRR13711360
SRR13711362
SRR13711363
可变剪切shell脚本
bash AS.sh
对于双端数据,只需要改变脚本中参考基因组和注释文件的下载地址和对应路径即可
#下载SRR文件
prefetch --option-file SRR_Acc_List.txt
mkdir project
cat SRR_Acc_List.txt |while read id;do mv ~/ncbi/public/sra/$id.sra project/ ;done #把下载的文件放到project里
cd project
wget https://ics.hutton.ac.uk/atRTD/RTD2/AtRTDv2_QUASI_19April2016.fa #下载参考基因组
wget https://ics.hutton.ac.uk/atRTD/RTD2/AtRTDv2_QUASI_19April2016.gtf #下载注释文件
gtf=/home/fanhz1/AS/homework/project/AtRTDv2_QUASI_19April2016.gtf #设置gtf路径
fa=/home/fanhz1/AS/homework/project/AtRTDv2_QUASI_19April2016.fa #设置fa路径
###上面的四行不同的物种需要自行更改
cd ../
#fastq-dump处理数据
mkdir fastq
cat SRR_Acc_List.txt | while read id ; do fastq-dump --split-3 --gzip project/$id.sra -O fastq ; done
#fastqc质控
mkdir qc
cat SRR_Acc_List.txt | while read id ; do fastqc -t 8 -o qc fastq/${id}_1.fastq.gz fastq/${id}_2.fastq.gz ; done
#使用trim_galore过滤,去接头
mkdir cleandata
cat SRR_Acc_List.txt | while read id ; do trim_galore --paired -q 25 --phred33 --length 36 --stringency 3 -o cleandata/ ${id}_1.fastq.gz ${id}_2.fastq.gz
#salmon建立参考基因组的索引
salmon index -k 25 -t $fa -i AtRTDv2.index
index=/home/fanhz1/AS/homework/AtRTDv2.index
#使用salmon计算tpm值
mkdir salmon
cat SRR_Acc_List.txt | while read id ; do salmon quant -i $index -l SF --gcBias --seqBias -1 cleandata/${id}_1_val_1.fq.gz -2 cleandata/${id}_2_val_1.fq.gz -p 8 -o salmon/${id}_output ; done
#使用salmon内置的脚本提取tpm值生成表达矩阵
mkdir tpm
multipleFieldSelection.py -i salmon/SRR*/quant.sf -k 1 -f 4 -o tpm/iso_tpm.txt
#简化样本id
cd tpm
perl -alne '{/(\|.*\|)\t/; ;s/$1//g;s/\|//g;print}' iso_tpm.txt | sed "s#_output##g"> iso_tpm_formatted.txt
cd ../
#对gtf文件使用suppa构建可变剪切事件综合文件
mkdir ioe
suppa.py generateEvents -i $gtf -o ioe/AtRTDv2.events -e SE SS MX RI FL -f ioe
#把所有可变剪接事件合并到一个文件夹里
cd ioe
awk 'FNR==1 && NR!=1 { while (/^<header>/) getline; } 1 {print}' *.ioe > AtRTDv2.all.events.ioe
cd ../
ioe=/home/fanhz1/AS/homework/ioe/AtRTDv2.all.events.ioe
#计算每个样本的psi值
mkdir psi
suppa.py psiPerEvent -i $ioe -e tpm/iso_tpm_formatted.txt -o psi/project_events
#对样本进行分组比较
cd psi
cut -f 1-5 project_events.psi > control.psi
cut -f 1-5 ../tpm/iso_tpm_formatted.txt > control.tpm
cut -f 1,6-9 project_events.psi > treat.psi
cut -f 1,6-9 ../tpm/iso_tpm_formatted.txt > treat.tpm
cd ../
#对两个分组的psi矩阵 和tpm矩阵进行 差异可变剪切事件寻找
suppa.py diffSplice -m empirical -gc -i $ioe --save_tpm_events -p psi/treat.psi psi/control.psi -e psi/treat.tpm psi/control.tpm -o project_diffSplice
网友评论