单端测序的RNAseq处理流程

srr_ftp_url=$1
srr_list=$2
my_current_env=$3
sra_dir=/media/whq/282A932A2A92F3D2/WHQ/mysradata/sra/srafile/sra
genome=/media/whq/282A932A2A92F3D2/WHQ/reference_genome/Homo_sapiens.GRCh38.dna.toplevel.fa
gtf=/media/whq/282A932A2A92F3D2/WHQ/reference_index/ensembl_GRCH38_hisat2_index/Homo_sapiens.GRCh38.99.gtf
index=/media/whq/282A932A2A92F3D2/WHQ/reference_index/ensembl_GRCH38_hisat2_index/ref.fa
get_count_trans=/media/whq/282A932A2A92F3D2/WHQ/myRscript/get_count_trans.R

#数据下载
#建环境，转移文件
#质控

ls -d SRR*|while  read id ;do input=$id/$id'.fastq.gz';output=$id/$id'.clean.fastq.gz';fastp -w 10 -i $input -o $output ;done

#比对、sort、删除sam

ls -d SRR*|while read id ;do input=$id/$id'.clean.fastq.gz';output=$id/$id'.sam';hisat2 -p 15 --dta -x $index -U $input  -S $output ;samtools sort -@ 10 --output-fmt BAM -o $id/$id'.sorted.bam' $id/$id'.sam' ;samtools index -@ 5 $id/$id'.sorted.bam';sam=$id/$id*'.sam';echo $sam;rm $sam ;done



#stringtie 这里希望能得到新的转录本！

ls -d SRR*|while read id;do input=$id/$id'.sorted.bam' ; output=$id/$id'.gtf' ;stringtie  -o $output -p 10 -G $gtf -B -l $id $input;echo $id ;done

ls -d  SRR*|while read id ;do find ./ -path './'$id'*' -name $id.gtf ;done >gtf.list

echo "merge function may error"
#注意，merge不知为何原因无法放在文件里跑，只能做代码跑

stringtie --merge  -p 10 -G $gtf  -o total_merged.gtf -l merge gtf.list

echo "merge function success!"

#使用gffcompare软件将merge与annotation比较
mkdir gffcompare
gff_prefix=./gffcompare/gffcompare
gffcompare -G -r $gtf -o $gff_prefix total_merged.gtf
mv gff*'map' gffcompare



#通过prepDE.py程序提取表达矩阵，stringtie2.1及以后就不会报错了，但是之前的数据要重新要stringtie跑一下
awk -F "/" '{print $2}' gtf.list >sample.list
ls -d  SRR*|while read id;do find ./ -path './'$id'*' -name *.merged.gtf ;done >merged.gtf1
paste sample.list merged.gtf1 >merged.gtf
#提取表达矩阵
mkdir count
mkdir transcript
prepDE.py -i merged.gtf -g ./count/gene_count_matrix_merged.csv -t ./transcript/transcript_count_matrix_merged.csv

#将产生的ball gown文件进行归类
mkdir ballgown_merged

ls -d SRR*|while read id ;do mkdir ./ballgown_merged/$id ;mv $id/*'ctab' ./ballgown_merged/$id/;done

#将数据读入R中用ballgown处理
#此处暂时不写！