这里是佳奥！经历了一个学期，做了几遍转录组，感受颇深。

原理没变，方法没变，但是思维需要转变，我们需要时刻保持清醒：

##这一步的代码是做什么

##我当前分析的是哪一个步骤

##这对于我的分析目的是否会有影响

做自己的项目不再是机械的重复步骤，而是要对自己的每一步负责，对课题组的经费支出负责。

那么，相信经过了先前的学习，我们已经掌握了转录组分析的全流程，这里，我会把使用的代码全部整理出来，分成上下篇，希望对做非模式物种的同学们有所帮助。

那么我们开始吧！

1 上传实验数据

##首先在服务器的目录下建立自己的目录
/mnt/disk/lja/project/

mkdir tree
cd tree

##我们后续所有的分析都在tree目录下

##上传数据推荐使用MobaXterm左侧文件栏的上传（向上箭头）功能
mkdir raw_fq
cd raw_fq

2 fastqc 质量控制

##首先激活conda小环境
conda activate rnaseq

##批量质控
ls *gz | xargs fastqc -t 10

##生成合成报告
multiqc ./

3 TrimGalore(依赖cutadapt) 过滤双端测序结果

##安装包下载，下载.gz至Linux下解压即可
https://github.com/FelixKrueger/TrimGalore/releases

##解压
tar -zxvf TrimGalore-0.6.7

##个人习惯，使用二进制版，且每次都要添加到环境变量
export PATH="$PATH:/mnt/disk/lja/biosoft/trimgalory/TrimGalore-0.6.7"

##批处理过滤

##先生成config文件
ls *R1.fastq.gz >1
ls *R2.fastq.gz >2
paste 1 2 > config

##设置文件所在路径
dir='/mnt/disk/lja/project/clean_fq'

##nohup挂起运行，可top查看进程情况
cat config | while read id
do
arr=($id)
fq1=${arr[0]}
fq2=${arr[1]}
echo $dir  $fq1 $fq2
nohup trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2 &
done

4 hisat2双端比对(这一步前备份质控后fq)

为什么这么说呢？因为比对的时候很吃内存和时间，中间要是有意外断了的话，fq文件容易损坏，就无法再继续分析了，就要回上一步重新过滤，为节省不必要的时间支出，尽量备份一遍。

让我们继续。

##到这一步，需要我们非模式种Tree的参考基因组文件，请找服务器管理员或者老师询问参考文件的位置

##建立索引，根据基因组大小，时间数小时不等
hisat2-build -p 8 Tree_V1.fasta genome

##原始文件名：Tree_1_val_1.fq.gz

##批量双端比对
cd /clean_fq
ls *.gz | cut -d "_" -f 1 | sort -u | while read id ; do
ls -lh ${id}_1_val_1.fq.gz  ${id}_2_val_2.fq.gz
nohup hisat2 -p 8 -x /mnt/disk/lja/refer/genome -1 ${id}_1_val_1.fq.gz -2 ${id}_2_val_2.fq.gz -S ${id}.hisat2.sam 
done

##批量转.bam
ls *.sam | while read id ; do (samtools sort -O bam -@ 5 -o $(basename ${id} ".sam").bam $id); done

##构建.bam索引（如果.bam太大会无法建立）
ls *.bam | xargs -i samtools index {}

5 subread (featureCounts)

##上有分析最后一步，统计count数，得到rawcount矩阵

##这一步需要我们非模式种Tree的.gtf注释文件

批量bam featureCounts:
gtf='/mnt/disk/lja/refer/Tree_V1.gtf'

nohup featureCounts -T 5 -p \
-a $gtf -o all.counts.txt \
/mnt/disk/lja/project/align/*.bam

至此，我们拿到了all.counts.txt，将它传输到我们的个人电脑，导入R开始下游分析吧！

6 写在后面

上游分析相对比较简单，虽然涉及的软件较多，但是思路是很明确的。

首先初步查看数据质量，进行过滤，再查看质量，开始与参考基因组比对，随后count比对数值。

上述代码适用于双端测序结果，以及有参考基因组的情况。无参情况请翻阅先前转录组的博客。

那么，上游分析有什么值得注意的呢？

1 文件规范命名

一个规范的命名可以节省很多不必要找文件的时间。尤其在Linux系统下，找文件需要更加多的操作。

##一个新项目一个目录
raw_fq：原始的fq文件
raw_qc：原始fq质控报告
clean_fq：过滤后的fq文件
clean_qc：过滤后的质控报告
align：比对（.sam .bam）
count：count矩阵

2 代码报错一步步排查

首先，文件是否存在，目录是否正确。

其次，conda环境是否正确激活，是否有该软件，该软件是否已被后台占用 。

最后，--help查看软件参数，是否是因为更新导致的参数变动。

3 固定流程

在没有出现重大分析问题的情况下，上游分析越快越好，因为我们更需要把时间花在下游分析上。

不仅是整理表达矩阵，还有在线注释蛋白.fa数据，以及不断筛选结果以对应我们的实验变量。

那么，我们非模式种转录组下游分析篇再见！