1：fastq检测测序质量，命令行如下
使用tmux将命令挂后台运行，不用担心因为网络问题掉线导致运行终止
tmux new -s fastqc
fastqc -q -t 3 -o /home/wanghongli/teacherW/BMK_DATA_20220325154115_1/Data/1/ *.fq.gz
-q：表示quiet，安静运行，最后生成报告
-t：表示线程，可以提高计算速度
-o：表示输出结果目录
运行结束会得到fastqc报告，
绿色表示通过检测，黄色警告，红色表示不通过，需要进一步处理原始reads。一般情况下，我们比较关注GC含量，Q20和Q30的比例以及是否存在接头（adaptor）、index以及其他物种序列的污染等。
2：去接头并质控
trim galore命令
tumx new -s trim_galore
trim_galore --illumina --fastqc --paired -o /home/wanghongli/teacherW/BMK_DATA_20220325154115_1/Data/1/ *.fastq.gz
--illumina：表示illumina接头
--fastqc：去完接头后进行fastqc质控
--paired：表示双端测序
-o：表示输出结果目录
输出结果：
生成的_val_fastq.gz文件是原始文件去完接头后得到的文件。通过网页打开质控文件，看看是否去接头成功。
3：比对
3.1.参考基因
小鼠 /home/wanghongli/mm10/mm10
人类 /home/wanghongli/hg38/hg38
3.2.基因注释
小鼠 /home/wanghongli/mm10/mm10.ncbiRefSeq.gtf
人类 /home/wanghongli/hg38/hg38.ncbiRefSeq.gtf
3.3. 比对 双端
小鼠 nohup hisat2 -p 16 -t -x /home/wanghongli/mm10/mm10 -1 T883-S01-01T0001_good_1.fq.gz -2 T883-S01-01T0001_good_2.fq.gz -S T883-S01-01T0001.sam &
人类 nohup hisat2 -p 16 -t -x /home/wanghongli/hg38/hg38 -1 T883-S01-01T0001_good_1.fq.gz -2 T883-S01-01T0001_good_2.fq.gz -S T883-S01-01T0001.sam &
3.4.sam转bam
nohup samtools view -bS T883-S01-01T0001.sam | samtools sort -@ 10 - -o T883-S01-01T0002.bam &
4：定量
4.1.构建bam索引
nohup samtools index T883-S01-01T0001.bam &
4.2.定量(两种定量方法)
4.2.1.
小鼠 featureCounts -T 10 -g gene_id -a /home/wanghongli/mm10/mm10.ncbiRefSeq.gtf -o T883-S01-01T0001_Counts.txt T883-S01-01T0001.bam
人类 featureCounts -T 10 -g gene_id -a /home/wanghongli/hg38/hg38.ncbiRefSeq.gtf -o T883-S01-01T0001_Counts.txt T883-S01-01T0001.bam
4.2.2.
小鼠 htseq-count -f bam -s no T883-S01-01T0002.bam /home/wanghongli/mm10/mm10.ncbiRefSeq.gtf > T883-S01-01T0002.count
人类 htseq-count -f bam -s no T883-S01-01T0002.bam /home/wanghongli/hg38/hg38.ncbiRefSeq.gtf > T883-S01-01T0002.count
5：无重复样本差异分析
5.1count转RPKM
RPKM=count/(length*文库大小）
每个样本中所有基因counts数求和是文库大小
5.2
5.2.1 安装gfold
Gfold软件是做没有生物学重复样本首选的软件
安装gsl最新版
以下所有路径均为解压后gsl内的路径
wget http://mirrors.ocf.berkeley.edu/gnu/gsl/gsl-latest.tar.gz
tar -zxv -f gsl-latest.tar.gz cd gsl-2.6/
./configure --prefix=/home/fanzhiyu_31946012_cyl_150/chayifenxi/gsl-2.6
make
make check(选做)
make install
安装gfold
cd .. wget https://zhanglab.tongji.edu.cn/softwares/GFOLD/gfold.V1.1.4.tar.gz
tar -zxf gfold.V1.1.4.tar.gz cd gfold.V1.1.4/
make
发现报错，输入以下命令
export CXXFLAGS="-g -O3 -I/home/fanzhiyu_31946012_cyl_150/chayifenxi/gsl-2.6/include -L/home/fanzhiyu_31946012_cyl_150/chayifenxi/gsl-2.6/lib"
export LD_LIBRARY_PATH="/home/fanzhiyu_31946012_cyl_150/chayifenxi/gsl-2.6/lib:"$LD_LIBRARY_PATH g++ -O3 -Wall -g main.cc -o gfold -lgsl -lgslcblas -I/home/fanzhiyu_31946012_cyl_150/chayifenxi/gsl-2.6/include -L/home/fanzhiyu_31946012_cyl_150/chayifenxi/gsl-2.6/lib
有warnning不用理会
./gfold -h
image
出现此安装成功
5.2.2.差异分析
5.2.2.1、把实验组和对照组的表达值分别建立两个csv文件，各包含如图五列，保证第一列顺序相同，如果不需要计算后两列的值，后两列可以随便填值
image
5.2.2.2、将自己生成的csv文件，再在 R 中统一一下格式（下面的代码是在简书作者没有猫但是有猫饼文章中粘的）
新生成的文件名我一般不改，就是 raw_data_Vec , raw_data_LSH 反正后面还会改掉
file1 = "ck.csv" raw_data_ck <- read.csv(file1, stringsAsFactors=FALSE)
file2 = "VI.csv" raw_data_vi <- read.csv(file2, stringsAsFactors=FALSE)
统一列名别出错
colnames(raw_data_ck) <- c("GeneSymbol", "GeneName", "Read Count", "Gene exon length", "RPKM") colnames(raw_data_vi) <- c("GeneSymbol", "GeneName", "Read Count", "Gene exon length", "RPKM")
导出文件，这个文件名我也不改，不然在 Linux 里改好麻烦
write.table(raw_data_ck, file="C:/Users/Administrator/Desktop/GFOLD/Sample1ck.read_cnt", row.names=F, col.names=F, quote=F, sep="\t")
write.table(raw_data_vi, file="C:/Users/Administrator/Desktop/GFOLD/Sample2vi.read_cnt", row.names=F, col.names=F, quote=F, sep="\")
5.2.2.3、将文件上传到服务器，运行如下命令
./gfold diff -s1 /home/fanzhiyu_31946012_cyl_150/chayifenxi/gfold.V1.1.4/Sample1ck -s2 /home/fanzhiyu_31946012_cyl_150/chayifenxi/gfold.V1.1.4/Sample2vi -suf .read_cnt -o sample1VSsample2.diff
5.2.2.4.筛选差异基因
大于0为上调，小于0下调（第3列是gfold值所在列，第4列是E值）
awk '{if($3>0 && $4=1) print $0}' sample1VSsample2.diff OFS='\t' > up_diff.gene
awk '{if($3<0 && $4=1) print $0}' sample1VSsample2.diff OFS='\t' > down_diff.gene
筛选两倍差异基因，即大于1和小于-1的
awk '{if($3>1 && $4=1) print $0}' sample1VSsample2.diff OFS='\t' > up_diff.gene-2
awk '{if($3<-1 && $4=1) print $0}' sample1VSsample2.diff OFS='\t' > down_diff.gene-2
网友评论