RNA-seq 分析流程（一）linux部分
基于miniconda3进行linux部分分析,由于很多软件是基于Python3.7写成的，建议使用Py3.7版本的miniconda3
miniconda3安装及注意事项见：
https://www.jianshu.com/p/914edc1de634
https://www.jianshu.com/p/42b7598fbc34
https://www.jianshu.com/p/e620115f7313
Step1:软件安装
数据资源下载，参考基因组及参考转录组(linux)
gtf ，genome，fa
质控，需要fastqc及multiqc(linux)
trimmomatic，cutadapt，trim-galore
对比（linux)
star,HISAT2,TOPHAT2, bowtie,bwa,subread
计数（linux）
featureCounts，htseq-counts,bedtools
nomalization 归一化，差异分析等（R 包）
DEseq2，edgeR，limma
注，fastqc和trim-galore需要下载openjdk，非常大，建议使用conda install --offline /home/xxx/src/openjdk-11.0.1-h516909a_1016.tar.bz2 线下安装
conda install -y fastqc trim-galore star bwa hisat2 bowtie2 subread tophat htseq bedtools deeptools salmon cutadapt multiqc sra-tools
conda install -c bioconda trimmomatic
Step2:原始数据校验
md5sum -c cdmd5.txt
Step3: 数据质控与矫正
3.1. fastqc
fastqc 批处理并将结果保存在qc文件夹下
fastqc -o ./qc/ *.fq.gz
multiqc整合fastqc结果
multiqc -o ./qc/ ./qc/*.zip
3.2. trim_galore 质控
去除接头
去除两端低质量碱基(-q 25)
最大允许错误率(默认-e 0.1)
去除<36的reads(--length 36)
切除index的overlap>3的碱基
reads去除以对为单位(--paired)
ls *_1.clean.fq.gz >1
ls *_2.clean.fq.gz >2
paste 1 2 >config
cat config
建立好config后，接下来就可以进行批量质控了
建立批量处理脚本
vim qc.sh
bin_trim_galore=trim_galore
mkdir filter-data
dir='/home/XXX/reads/filter-data'
cat config |while read id
do
arr=($id)
fq1=${arr[0]}
fq2=${arr[1]}
nohup $bin_trim_galore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2 &
done
拿到config和qc.sh后
bash qc.sh config
注意trim_galore是平行批量处理，对于在自己电脑上做比对的人来说，对电脑伤害比较大
3.3. 异常碱基剪切
发现本次测序数据的前15个碱基碱基比例异常，因此决定剪掉起始的25bp碱基
vim trim.sh
ls *.gz
ls *.gz|while read id;do echo $id;
cutadapt -u 25 -o /home/XXX/trim.$id /home/XXX/$id;
done
去掉样本名称前的trim
rename "s/trim.//" *fq.gz
Step4: hisat2 比对 (可选1）
hisat2 build index
hisat2_extract_exons.py tair10_ensemble_chr.gtf >exon_arab.txt
hisat2_extract_splice_sites.py tair10_ensemble_chr.gtf >ss_arab.txt
hisat2-build -p 2 --exon /home/polya/Public/genome/Arab/tair/exon_arab.txt --ss /home/polya/Public/genome/Arab/tair/ss_arab.txt /home/XXX/Public/genome/Arab/tair/tair10_ensemble_chr.fa /home/XXX/Public/genome/Arab/index/hisat/index
hisat2 mapping 使用vim map.sh, 构建比对，可以一个一个比对，也可以写个循环比对
hisat2 -p 2 -x /home/polya/Public/genome/Arab/index/hisat/index -1 /home/polya/data/clean/trim/sample_read1.clean.fq.gz -2 /home/polya/data/clean/trim/sample_read2.clean.fq.gz -S /home/polya/data/clean/sam/sample.hisat.sam;
Step4: STAR 比对 (可选2）
star build the index
STAR --runThreadN 6 --runMode genomeGenerate --genomeDir /home/polya/Public/genome/Arab/tair10/ --genomeFastaFiles /home/polya/Public/genome/Arab/tair10/TAIR10_chr_all.fas --sjdbGTFfile /home/polya/Public/genome/Arab/tair10/TAIR10_GFF3_genes.gff --sjdbGTFfeatureExon exon --sjdbGTFtagExonParentTranscript ID --sjdbGTFtagExonParentGene Parent
使用vim map.sh, 构建比对，可以一个一个比对，也可以写个循环比对
STAR --runThreadN 12 --readFilesCommand zcat --alignEndsType EndToEnd --outSAMtype SAM --outFilterMultimapNmax 1 --genomeDir /home/polya/Public/genome/Arab/tair10/ --readFilesIn rep1.R1.clean_val_1.fq.gz rep1.R2.clean_val_2.fq.gz --outFileNamePrefix /home/polya/Public/data/sam/rep1;
Step5: sam to bam, 继续使用vim xxx.sh,使用chmod 777 给sh权限，然后nohup ./xxx.sh 2>&1&后台挂起，下同
ls *.sam
ls *.sam|while read id;do echo $id;
samtools view -Sb $id > ${id%%.*}.bam;
done
Step6:bam to sorted bam
ls *.bam
ls *.bam|while read id;do echo $id;
samtools sort $id -o ${id%%.*}.sorted.bam;
done
rename 's/Aligned.out//g' *
Step7: generate flagstat for summary of mapping
ls *sorted.bam |while read id ;do ( samtools flagstat -@ 1 $id > $(basename ${id} ".bam").flagstat );done
index for IGV visualization
ls *sorted.bam|xargs -i samtools index {} #构建索引，拿到IGV里面看
multiqc ./*.flagstat
Step8 featurecounts to calculate gene expression,常规RNA-seq分析用这个就行 （多数情况选择）
ls *sorted.bam
ls *sorted.bam|while read id;do echo $id;
featureCounts -T 6 -p -t exon -g gene_id -a /home/polya/Public/genome/Arab/tair10/TAIR10_ensemble_Chr.gtf -o all.gene.txt *.bam;
done
这一步结束就可以移动到R 中进行下游的差异分析了：详见RNAseq分析流程（二）https://www.jianshu.com/p/7c0d61363ce3
少数情况：
Step8: #featurecounts to calculate first exon expression，针对RNA 加工中的转录本差异用的，小众 （仅少见情况选择）
ls *sorted.bam
ls *sorted.bam|while read id;do echo $id;
featureCounts -T 6 -t exon -g Parent -a /home/polya/Public/genome/Arab/tair10/TAIR10_GFF3_fexon.gff -o all.exon.txt *.bam;
done
如果想使用bedgraph看IGV
Step9： bam to bedgraph, Note: plus is actual minus for strand-specific reads，这部分是生成IGV的bedgraph文件
ls *.sorted.bam
ls *.sorted.bam|while read id;do echo $id;
genomeCoverageBed -ibam $id -bga -split -g /home/polya/Public/genome/Arab/tair10/TAIR10_chr_all.fas > ${id%%.*}.bedgraph
done
