基于miniconda3進(jìn)行l(wèi)inux部分分析,由于很多軟件是基于Python3.7寫成的，建議使用Py3.7版本的miniconda3
miniconda3安裝及注意事項(xiàng)見：
http://www.itdecent.cn/p/914edc1de634
http://www.itdecent.cn/p/42b7598fbc34
http://www.itdecent.cn/p/e620115f7313
Step1:軟件安裝
數(shù)據(jù)資源下載，參考基因組及參考轉(zhuǎn)錄組(linux)
gtf ，genome，fa
質(zhì)控，需要fastqc及multiqc(linux)
trimmomatic，cutadapt，trim-galore
對比（linux)
star,HISAT2,TOPHAT2, bowtie,bwa,subread
計(jì)數(shù)（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:原始數(shù)據(jù)校驗(yàn)
md5sum -c cdmd5.txt
Step3: 數(shù)據(jù)質(zhì)控與矯正
3.1. fastqc
fastqc 批處理并將結(jié)果保存在qc文件夾下
fastqc -o ./qc/ *.fq.gz
multiqc整合fastqc結(jié)果
multiqc -o ./qc/ ./qc/*.zip
3.2. trim_galore 質(zhì)控
去除接頭
去除兩端低質(zhì)量堿基(-q 25)
最大允許錯誤率(默認(rèn)-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后，接下來就可以進(jìn)行批量質(zhì)控了
建立批量處理腳本
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. 異常堿基剪切
發(fā)現(xiàn)本次測序數(shù)據(jù)的前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, 構(gòu)建比對，可以一個一個比對，也可以寫個循環(huán)比對
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, 構(gòu)建比對，可以一個一個比對，也可以寫個循環(huán)比對
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, 繼續(xù)使用vim xxx.sh,使用chmod 777 給sh權(quán)限，然后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 {} #構(gòu)建索引，拿到IGV里面看
multiqc ./*.flagstat
Step8 featurecounts to calculate gene expression,常規(guī)RNA-seq分析用這個就行 （多數(shù)情況選擇）
ls *sorted.bam
ls *sorted.bam|while read id;do echo $id;
featureCounts -T 6 -p -t gene -g gene_id -a /home/polya/Public/genome/Arab/tair10/TAIR10_ensemble_Chr.gtf -o all.gene.txt *.bam;
done
這一步結(jié)束就可以移動到R 中進(jìn)行下游的差異分析了：詳見RNAseq分析流程（二）http://www.itdecent.cn/p/7c0d61363ce3
少數(shù)情況：
Step8: #featurecounts to calculate first exon expression，針對RNA 加工中的轉(zhuǎn)錄本差異用的，小眾 （僅少見情況選擇）
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
