上游分析的最后一部分，这里是佳奥，让我们继续吧！

1 获取实验数据（双端测序）

raw_fq：下载sra文件，转fq文件。

clean_fq：trim_galory

ls *_1.fastqc.gz >1
ls *_2.fastqc.gz >2
paste 1 2 > config

dir='../clean_fq'
cat $config_file | while read id
do
arr=($id)
fq1=${arr[0]}
fq2=${arr[1]}
echo $dir  $fq1 $fq2
trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2 &
done

qc：质量控制，并查看.html报告。

2 hisat2比对

2.1 建立索引

##网上下载或者自行构建

##构建索引
hisat2-build -p 16 Arabidopsis_thaliana.TAIR10.28.dna.genome.fa genome

align：比对后的.bam/.sam文件。

2.2 批量比对

index=../hisat2index
fq1=..
fq2=..

sample=${arr[0]}
if((i%$number1==$number2))
    then
        if[! -f $sample.bam]; then
            start=$(date +%s.%N)
            echo hisat2 'date'
            hisat2 -P 4 -x $index -1 $fq1 -2 $fq2 | samtools sort -@ 4 -o $sample.bam
            echo hisat2 'date'
            dur=$(echo "$(date +%s.%N) - $start" | bc)
            printf "Execution time for hisat2 : %.6f seconds" $dur  
        fi##endforffiles
    fi##end for number1
    i=$((i+1))
done ##end for $config- file

3 bamCoverage转bw文件

##方法一
analysis_dir=$1
config_file=$2
number1=$3
number2=$4

cat $config_file | while read id
do
echo $id
    echo $id
    file=$(basename $id )
    sample=${file%%.*}
    echo $sample
    
    if((i%$number1==$number2))
    then 
        if [ ! -f $sample.bw ]; then
            start=$(date +%s.%N)
            echo bamCoverage 'date'
            bamCoverage -b $id -o $sample.bw --normalizeUsing RPKM -p 4
            echo bamCoverage 'date'
            dur=$(echo "$(date +%s.%N) - $start" | bc)
            printf "Execution time for bamCoverage : %. 6f seconds" $dur
        fi##end for if files
    fi##end for number1
    i=$((i+1))
done##end for $config_file

##方法二、一次性全部提交，如果是公共服务器等着管理员来联系你（所以会有方法一的方法）
ls *bam| while read id; do(bamCoverage -b $id -o ${id/.merge.bam/.bw} --normalizeUsing RPKM -p 4 &);done

4 featureCounts差异分析

##方法一、批量bam featureCounts
gtf='/home/kaoku/rnaseq/biotree_plant/refer/Arabidopsis_thaliana.TAIR10.28.gtf.gz'

featureCounts -T 5 -p \
-a $gtf -o all.counts.txt \
/home/kaoku/rnaseq/biotree_plant/data/sam_bam_bai/*.bam

##方法二
featureCounts -T 4 -p -t exon -g gene_name -a $gtf -o all.counts.id.txt ../bams/*.bam 1>counts.id.log 2>&1

最后会获得一个矩阵文件，all.id.counts.txt，可以再multiqc ./看一下结果。

至此上游分析结束。

后续的内容便是把表观调控的图表重现出来。

我们下一篇再见！