🍊码农

HISAT2_Stringtie_DEseq2 shell流程

2018-05-22  本文已影响102人  547可是贼帅的547

#!/bin/sh

#参数的输入

read  -a control -p "请输入control组文件名,空格分开:"

echo -e "\n"

read -a treat -p "请输入treatment组文件名,空格分开:"

echo -e "\n"

read  -p "请输入要使用的线程数:" thread

read  -p "请输入参考基因组的索引:" index

echo -e "\n"

read  -p "请输入参考基因组gtf文件名:" gtf

echo -e "\n"

read  -p "请输入R脚本文件名:" Rscript

echo -e "\n"

echo "control组文件:  ${control[*]} "

echo "treatment组文件:  ${treat[*]}  "

echo "线程数: $thread"

echo "参考基因组的索引: $index"

echo "参考基因组gtf文件: $gtf"

#echo "参考基因组fasta文件: $fasta"

echo -e "\n"

#比对部分

echo "control组mapping、alignment"

echo -e "\n"

rm mergelist.txt

for ((i=0;i<${#control[*]};i++));do

echo ${control[$i]}

echo -e "\n"

hisat2 --dta -x $index -p $thread -U ${control[$i]} -S ${control[$i]}.sam

samtools sort -@ $thread -o ${control[$i]}.bam ${control[$i]}.sam

samtools index -@ $thread ${control[$i]}.bam

stringtie -p $thread -G $gtf -o ${control[$i]}.gtf -A ${control[$i]}.tab -B -e -l ${control[$i]} ${control[$i]}.bam

echo -e "${control[$i]}\t./${control[$i]}.gtf" >> mergelist.txt

done

echo "treatment组mapping、alignment"

echo -e "\n"

for ((i=0;i<${#treat[*]};i++));do

echo ${treat[$i]}

echo -e "\n"

hisat2 --dta -x $index -p $thread -U ${treat[$i]} -S ${treat[$i]}.sam

samtools sort -@ $thread -o ${treat[$i]}.bam ${treat[$i]}.sam

samtools index -@ $thread ${treat[$i]}.bam

stringtie -p $thread -G $gtf -o ${treat[$i]}.gtf -A ${treat[$i]}.tab -B -e -l ${treat[$i]} ${treat[$i]}.bam

echo -e "${treat[$i]}\t./${treat[$i]}.gtf" >> mergelist.txt

done

echo "转换为DEseq2可用的matrix"

python ~/Software/prepDE.py -i mergelist.txt

echo "运行指定的R脚本运行DEseq2"

Rscript $Rscript ${#control[*]} ${#treat[*]}

上一篇下一篇

猜你喜欢

热点阅读