上游分析的最后一部分，這里是佳奧，讓我們繼續(xù)吧！

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 ./看一下結果。

至此上游分析結束。

后續(xù)的內容便是把表觀調控的圖表重現出來。

我們下一篇再見！