从下载数据开始完整跑hisat2-stringtie-提取表达矩

参考文章：https://www.jianshu.com/p/1f5d13cc47f8
https://davetang.org/muse/2017/10/25/getting-started-hisat-stringtie-ballgown/

第一部分代码：

#!/bin/bash
#人-转录组-hisat2+stringtie
#参数准备（注意，设置为自己账号下的路径）
#注意，输入的参数全部为绝对路径！！！
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

#从ENA中下载fastq序列
cd $sra_dir
rm ascp_com2
cat $srr_ftp_url|while read id;do id=${id#*vol1/} ; echo "ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/${id} ./">>ascp_com2 ; done
echo "已经建立好下载通道"
./ascp_com2


#将下载好的fastq数据复制到my_current_env
cd $my_current_env &&mkdir fastq&&cd fastq
cd $sra_dir
cat $srr_list |while read id ;do read1=`ls $id'_1'*` ;read2=`ls $id'_2'*` ; cp $read1 $my_current_env/fastq/ ;cp $read2 $my_current_env/fastq/ ;done
cd -
#fastq中新建文件夹并把fastq文件移动进去
ls SRR*_1*|while read id ;do id=${id%%_1*};mkdir $id ;file1=`ls $id"_1"*`;file2=`ls $id"_2"*`;mv $file1 $id/;mv $file2 $id/;done
#使用fastp质控

echo "开始质控"
ls -d SRR*|while read id ;do read1=`ls $id/$id'_1'*` ; read2=`ls $id/$id'_2'*`  ;output1=$id/$id'_1.clean.fastq.gz' ; output2=$id/$id'_2.clean.fastq.gz'  ;echo $read1 $output1 $read2 $output2 ;fastp -i $read1 -o $output1 -I $read2 -O $output2 ;done
#使用fastqc进行质控
mkdir fastqc
ls -d SRR*|while read id ;do read1=$id/$id'_1.clean.fastq.gz';read2=$id/$id'_2.clean.fastq.gz';fastqc -o fastqc $read1 $read2;done


#比对、sort、删掉bam


ls -d SRR*|while read id ;do read1=$id/$id'_1.clean.fastq.gz' ; read2=$id/$id'_2.clean.fastq.gz' ;hisat2 -p 15 --dta -x $index -1 $read1 -2 $read2 -S $id/$id'.sam' ;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处理
#此处暂时不写！

第三部分代码：只想对annotation.gtf中已知的转录本进行定量

#stringtie 这里只想得到参考的gtf中的转录本定量！

ls -d SRR*|while read id;do stringtie -e -A $id/$id'_gene_abund.tab' -C $id/$id'_cov_refs.gtf' -B -p 10 -G $gtf -o $id/$id'.annotated.gtf' $id/$id'.sorted.bam' ;done

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

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

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

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

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