这里是佳奥！经历了一个学期，做了几遍转录组，感受颇深。
原理没变，方法没变，但是思维需要转变，我们需要时刻保持清醒：
##这一步的代码是做什么
##我当前分析的是哪一个步骤
##这对于我的分析目的是否会有影响
做自己的项目不再是机械的重复步骤，而是要对自己的每一步负责，对课题组的经费支出负责。
那么，相信经过了先前的学习，我们已经掌握了转录组分析的全流程，这里，我会把使用的代码全部整理出来，分成上下篇，希望对做非模式物种的同学们有所帮助。
那么我们开始吧！
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数据，以及不断筛选结果以对应我们的实验变量。
那么，我们非模式种转录组下游分析篇再见！
网友评论