Transcription analysis 183个杜洛克猪的rna-seq数据
在linux里使用后台命令进行下载，虽然上传的wget文件里都有，但是下载的少了两个，利用R语言的setdiff函数比对原始数据和下载的数据就可以很快找到缺失的两个数据（match.R）。
nohup wget -O "NAME" -c site>NAME.out & #download
在下载数据的时候发现，NCBI上有的是有fastq格式的数据的有的没有，所以统一下载NCBI格式的数据。 首先将数据名集中到一个文件夹里，按每十行来拆分成若干name文件，使用集群推荐的脚本进行处理。
ls S* E* >test #汇集S、E开头的文件到test文件夹里
split -| 10 test -d -a2 name_ #拆分
script.sh #集群脚本
#!/bin/bash
#SBATCH --job-name=RNA ##RNA
#SBATCH --partition=big ##作业申请的分区名称
#SBATCH --nodes=1 ##作业申请的节点数
#SBATCH --ntasks-per-node=1 ##作业申请的每个节点使用的核心数
#SBATCH --error=chaifen.err
#SBATCH --output=chaifen.out
echo "process will start at : "
date
echo "++++++++++++++++++++++++++++++++++++++++"
cd /public/agis/liuyuwen_group/wangchao/double
conda init bash
conda activate RNA-seq
cat name_00 |while read id;do fastq-dump --split-3 $id;done # 拆分数据
date
sbatch scrpit.sh #提交任务
还有一步很重要是通过NCBI的数据来判断下载的数据是否是链特异性建库，因为链特异性建库可以区分转录本来自正义反义链的，如果不区分的化直接QC然后比对可能会造成数据的缺失。 接下来进行数据的质控，主要是去除adapter和N碱基多的reads去除低质量的片段，最后保留duplication进行下一步的定量。 在集群里面提交脚本
#!/bin/bash
#SBATCH --job-name=RNA ##RNA
#SBATCH --partition=low ##作业申请的分区名称
#SBATCH --nodes=1 ##作业申请的节点数
#SBATCH --ntasks-per-node=1 ##作业申请的每个节点使用的核心数
#SBATCH --error=sample_0.err
#SBATCH --output=sample_0.out
echo "process will start at : "
date
echo "++++++++++++++++++++++++++++++++++++++++"
cd /public/agis/liuyuwen_group/wangchao/double/unspecificity
cat sample__0 |while read id;
do echo $id
arr=($id)
fq2=${arr[2]}
fq1=${arr[1]}
sample=${arr[0]}
trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired --gzip -o clean_data $fq1 $fq2
done
date
后面就会生成质控之后的文件
下面详细记录一下脚本里面的语法
--quality：设定Phred quality score阈值，默认为20。
--phred33：：选择-phred33或者-phred64，表示测序平台使用的Phred quality score。
--adapter：输入adapter序列。也可以不输入，Trim Galore!会自动寻找可能性最高的平台对应的adapter。自动搜选的平台三个，也直接显式输入这三种平台，即--illumina、--nextera和--small_rna。
--stringency：设定可以忍受的前后adapter重叠的碱基数，默认为1（非常苛刻）。可以适度放宽，因为后一个adapter几乎不可能被测序仪读到。
--length：设定输出reads长度阈值，小于设定值会被抛弃。
--paired：对于双端测序结果，一对reads中，如果有一个被剔除，那么另一个会被同样抛弃，而不管是否达到标准。
--retain_unpaired：对于双端测序结果，一对reads中，如果一个read达到标准，但是对应的另一个要被抛弃，达到标准的read会被单独保存为一个文件。
--gzip和--dont_gzip：清洗后的数据zip打包或者不打包。
--output_dir：输入目录。需要提前建立目录，否则运行会报错。
-- trim-n : 移除read一端的reads
trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired --gzip -o clean_data $fq1 $fq2 #所以这个脚本大致是使用phred进行质量微调，丢弃长度为35bp的reads，与修剪序列所需的adapter序列重叠、重叠数为4，错误率小于0.1，--gzip：清洗后的数据zip打包。
首先下载参考基因组 野猪的基因组信息 然后需要对于之前质控的文件，进行参考基因组的比对。 提交脚本
#!/bin/bash
#SBATCH --job-name=RNA ##RNA
#SBATCH --partition=low ##作业申请的分区名称
#SBATCH --nodes=1 ##作业申请的节点数
#SBATCH --ntasks-per-node=8 ##作业申请的每个节点使用的核心数
#SBATCH --error=sample_0.err
#SBATCH --output=sample_0.out
echo "process will start at : "
date
echo "++++++++++++++++++++++++++++++++++++++++"
cd /public/agis/liuyuwen_group/wangchao/double/unspecificity/clean_data #这里是非特异性建库的质控后的数据的地址
reference=/public/agis/liuyuwen_group/wangchao/pig_refrence/Sus #参考基因组的地址
cat sample__0 |while read id;
do echo $id
arr=($id)
fq2=${arr[2]}
fq1=${arr[1]}
sample=${arr[0]}
hisat2 -x $reference -q -1 $fq1 -2 $fq2 --new-summary | samtools sort -O bam -@ 5 -o - > BAM/${sample}.bam
#**如果是链特异性需要加上参数--rna-strandness RF --rna-strandness FR its paired end data and forward stranded，设置为FR，反之为RF
done
date
比对软件有tophat，hisat2，STAR，这里记录一下hisat2的参数
# **注意连特异性建库的时候要加上--rna-strandness RF这个参数**
这里使用的featucount软件，RNA-seq的计数较为复杂，因为需要考虑到外显子剪切。一种计数方法是数一下与每一个被注释的外显子重合的read，另外一种方法是数一下与每一个基因区域重合的read。 提交脚本
#!/bin/bash
#SBATCH --job-name=RNA ##RNA
#SBATCH --partition=low ##作业申请的分区名称
#SBATCH --nodes=3 ##作业申请的节点数
#SBATCH --ntasks-per-node=8 ##作业申请的每个节点使用的核心数
#SBATCH --error=sample_1.err
#SBATCH --output=sample_1.out
echo "process will start at : "
date
echo "++++++++++++++++++++++++++++++++++++++++"
cd /public/agis/liuyuwen_group/wangchao/double/unspecificity/clean_data/BAM
featureCounts -p -g gene_id -t exon -a Sus_scrofa.Sscrofa11.1.103.gtf -o featureCounts_results *.bam #p代表双端， -t表示feature-type选择外显子，-a表示参考gtf文件名，支持Gzipped文件格式，-g：当参考的gtf提供的时候，我们需要提供一个id identifier 来将feature水平的统计汇总为meta-feature水平的统计，默认为gene_id，注意！选择gtf中提供的id identifier
date
# 注释文件添加绝对路径， 特意连比对加 -s
# 将用到的每个参数的含义写到下方
# 所有成功运行的 .err .out文件不要删除，后面需要整合一些结果
#-p 表示针对双端测序数据 ，-g指定统计的meta-feature一般是gene_id，-t指定统计的feature一般是外显子，-o输出为，-a 输出GTF\GFF的注释文件，输入前切记要添加绝对路径。
使用R软件对featureCount 的结果进行处理
#! /usr/bin/env Rscript
## Tpm
tpm <- function(counts, lengths) {
rate <- counts / lengths
rate / sum(rate) * 1e6
}
## read table from featureCounts output
args <- commandArgs(T) ##commandArg是R里面传递参数的函数如arg[1],arg[2]就是传递脚本里面的第一个参数，脚本里面有一个txt
tag <- tools::file_path_sans_ext(args[1]) ##tools::file_path_sans_ext：r内置函数该文件包可获取不带扩展名的文件
ftr.cnt <- read.table(args[1], sep="\t", stringsAsFactors=FALSE, ##读取表
header=TRUE)
library(dplyr)
library(tidyr)
ftr.tpm <- ftr.cnt %>% ##%>%是管道操作，来自dplyr包的管道函数，其作用是将前一步的结果直接传参给下一步的函数，从而省略了中间的赋值步骤
gather(sample, cnt, 7:ncol(ftr.cnt)) %>% 查看featureCount第七列是SRA开头的文件，，gather函数是宽转长的函数，这里是7至后面的列作为cnt，其他的作为sample
group_by(sample) %>% group_by是以XX为标准排列的函数
mutate(tpm=tpm(cnt, Length)) %>% mutate是对已有列为运算，添加为新列
select(-cnt) %>% 选择除cnt之外的列
spread(sample, tpm) # 长转宽
write.table(ftr.tpm, file=paste0(tag, "_TPM.txt"), sep="\t", row.names=FALSE, quote=FALSE)
在脚本中提交这个Rscript，完成定量。