RNA-seq全细胞操作流程

# 首先创建数据文件夹
mkdir RNA_seq_all
mkdir origin_data
# 在origin_data下下载原始数据

# 由于这次数据下载量比较大，我尝试写一个脚本自动下载
# 需要下载的数据集在这个网址
https://www.ncbi.nlm.nih.gov//geo/query/acc.cgi?acc=GSE72379

# 开始就遇到了我的第一个问题，在重复上一个流程的时候，每个状态的细胞首先分为两个重复组，然后在每个重复组里又有两个重复，相当于每个形态的细胞有四组RNA-seq数据，是选择一份还是全部都下载呢，这里我想先每个细胞只拿两组数据跑一下
# 另外在Germinal Vesicle时期中只有两组细胞，所以这里我们全都用上

# 开始脚本

# 首先写了一个ID的list
# 这个应该...不用教吧，搜一下vim
SRR2183357 
SRR2183358
SRR2183359
SRR2183361
SRR2183363
SRR2183365
SRR2183367
SRR2183369
SRR2183371
SRR2183373
SRR2183375
SRR2183377
SRR2183379
SRR2183381
SRR2183383
SRR2183385


# 测试一下能否逐行读取数据，注意这里是``而不是''，在键盘esc的下面那个按键，反向引号
for ID in `cat SRR_ID`;do echo $ID; done

# 成了，接下来运行看看
for ID in `cat SRR_ID`; do nohup fasterq-dump --progress $ID & done

# 刚开始是成功的，但是后面就全都崩了
# 错误代码
fasterq-dump quit with error code 3

# 只能和上次一样一个一个下载了，就是比较浪费时间
# 或者在写shell的时候考虑完成了当前的下载再执行下一个?
for ID in `cat SRR_ID`; do fasterq-dump --progress $ID ; done
# 成功了，要留600G左右的硬盘空间才可以

# 查看当前目录所占大小
du -h --max-depth=1

# 或者
ls -lht

# 下一步尝试给文件批量命名，并提取出前100万行数据进行预分析顺便测试脚本的可行性，之后换原始数据进行全局比对
# 创建测试目录
mkdir test_data

# 说来惭愧，这部分我就想了半个上午，换了几种不同的方案，最后还是得自己手动写一个文件名目录再半自动批量命名，还得多多学习
# 文件名my_name
GV1_1 GV1_2 GV2_1 GV2_2 MI1_1 MI1_2 MI2_1 MI2_2 MII1_1 MII1_2 MII2_1 MII2_2 PN1_1 PN1_2 PN2_1 PN2_2 CL1_1 CL1_2 CL2_1 CL2_2 MO1_1 MO1_2 MO2_1 MO2_2 IC1_1 IC1_2 IC2_1 IC2_2 TP1_1 TP1_2 TP2_1 TP2_2

# 脚本代码
i=1;for file in *.fastq;do name=$(cat my_name|awk '{printf $"'$i'"}'); cat $file|head -4000000 > ../test_data/${name}.fastq; i=`expr $i + 1`; done

# 脚本讲解
# 本来以为这里会很简单，但没想到在不同的地方卡了很久，首先是最简单的变量赋值，一定注意在shell脚本中的赋值，在运算符号两边需要空格，且需要`expr`，注意并不是引号是反引号，此外还有在for中嵌套正则表达式cat|awk，尤其是关于awk变量的使用方法在$i还要再加个$,但这里还有个小bug，就是在编辑文件目录的时候，需要按照数字的顺序写文件名，否则顺序会对错，我在运行之后抽测了几个数据发现没什么问题，但是这里不注意会是一个隐患，之前考虑到使用if或者switch去做切换，但这样代码就写的又臭又长，太不美观了
# 最后这个代码会输出在根目录下的test_data文件夹，注意地址

# -t调用两个核 -o输出到哪个文件 -q沉默模式
fastqc -t 2 -o ../FastQC_result/1_fastqc/ -q *.fastq &
# 注意输出目录，需要自己新建

# 如果是服务器，可以使用filezilla把数据文件down下来查看

# -q 20 去除质量分数低于20的碱基序列，另外去除--nextera的引物
trim_galore -q 20 --nextera --paired *.fastq  -o ../trim_report/ &

# 这里把txt报告放到一个文件夹里方便批量处理文件
mkdir report
mv *.txt report

# 写脚本把这个软件写的名字给改成自己喜欢的
for file in *.fq;do name=$(ls $file|cut -d_ -f1,2);mv $file trim_${name}.fq;done

# 这样就改成了类似trip_CL1_1.fq的格式了，方便知道这是做了哪一步的文件
# 这里想用fastqc再看看是否处理了引物，还有其他地方有没有什么不同
# 注意输出文件夹
fastqc -t 2 -o ../FastQC_result/20220117_2_trim_fastqc/ -q *.fq &

# 这里rRNAdust的作用是除去冗杂的rRNA,貌似这个测序Helicos数据集的错误率很高，所以要把错误丢弃
# 安装rRNAdust
wget https://fantom.gsc.riken.jp/5/suppl/rRNAdust/rRNAdust1.06.tgz
tar -zxvf rRNAdust1.06.tgz
rm rRNAdust1.06.tgz

# cd到解压文件夹目录
cd rRNAfilter
make

# 这里用实验室服务器的话会遇到一个bug，提示无权限，你只要在自己的文件夹下新建一个bin然后导入PATH就好了
mkdir $HOME/bin
export PATH=$PATH:$HOME/bin
cp rRNAdust $HOME/bin

# 使用脚本进行批量处理
for file in *.fq;do name=$(ls $file|cut -d_ -f2,3);cat ${file} | rRNAdust -t 8 ../../U13369.1.fa > ../rRNAdust_result/rRdust_${name};done

# 再查看一次fastqc
fastqc -t 2 -o ../FastQC_result/20220117_3_rRNAdust_fastqc/ -q *.fq &

# 使用seqkit比较两个fastq文件，这里因为使用rRNAdust把污染的部分去除了，因为这个数据是双端测序，所以去除了一端，另一端也要去除
for file in *.fq;do name=$(ls $file|cut -d_ -f1,2);seqkit pair -1 ${name}_1.fq -2 ${name}_2.fq;mv *.paired.fq ../seqkit_result/;done

# 让我纳闷的是seqkit不能指定输出文件夹吗，也可能是我没看到，不过就脚本加一行的事情，顺便我们再改个名字
for file in *.fq;do name=$(ls $file|cut -d_ -f2,3);mv ${file} ${name};done

# 在服务器上运行mapping

qsub qsub.STAR.sh
qsub -I

#!/bin/bash
#PBS -q MEDIUM
#PBS -l select=1:ncpus=8
#PBS -V
cd ~/DATA/RNA_seq_1/test_data/seqkit_result/

fa="/user2/tanpaku/mizuguchi/HuYajie/DATA/RNA_seq_1/genomeDIr/hg38_index_dir/"
cpu=8
for file in *.fq;do
        name=$(ls $file|cut -d_ -f1);
        STAR --readFilesIn ${name}_1.paired.fq ${name}_2.paired.fq --runThreadN $cpu --genomeDir $fa --outSAMunmapped Within --outSAMtype BAM SortedByCoordinate --outFileNamePrefix ../STAR_result/${name};
done;

# 由于第一次上传到服务器运行数据，所以有很多不懂的地方，注意这里的--genomeDir地址需要输入绝对地址，输入相对地址可能找不到

# 使用featureCount计算表达量
# 这一步可以把所有数据生成到一个txt文件中
# 注意文件路径
featureCounts -T 6 -p -t exon -g gene_id -a ~/DATA/RNA_seq_1/genomeDIr/gencode.v38.annotation.gtf -o Counts.txt *.bam

library(edgeR)
mydata <- read.table("C:/Lab/DATA/20220118_Counts.txt",header = TRUE, quote = '\t', skip = 1)
sampleNames <- c("CL1","CL2","GV1","GV2","IC1","IC2","MI1","MI2","MII1","MII2","MO1","MO2","PN1","PN2","TP1","TP2")
names(mydata)[7:22] <- sampleNames
countMatrix <- as.matrix(mydata[7:22])
group <- factor(c("CL","CL","GV","GV","IC","IC","MI","MI","MII","MII","MO","MO","PN","PN","TP","TP"))
y <- DGEList(counts = countMatrix,group = group)
keep <- rowSums(cpm(y)>1) >= 2
y <- y[keep,,keep.lib.sizes=FALSE]
y <- calcNormFactors(y)
design <- model.matrix(~group,y)
install.packages('statmod')
y <- estimateDisp(y,design,robust = TRUE)
y$common.dispersion
plotBCV(y)
fit = glmFit(y, design)
lrt = glmLRT(fit,coef=2)
results = topTags(lrt, adjust.method="fdr", n=dim(y)[1])$table
key = rownames(results)
logcpm = cpm(y, log=TRUE)
results = cbind(logcpm[key,], results$logFC, results$logCPM, FDR=results$FDR)
result <- results
result <- data.frame(result)
library(ggplot2)
pairs(~GV1+GV2+MI1+MI2+MII1+MII2+PN1+PN2+CL1+CL2+MO1+MO2+IC1+IC2+TP1+TP2,data = result)
# 在进行最后一步的时候需要时间，如果R卡死了你可以试着减少变量操作

# download the raw data
for ID in `cat SRR_ID`; do fasterq-dump --progress $ID ; done

# filename:my_name
GV1_1 GV1_2 GV2_1 GV2_2 MI1_1 MI1_2 MI2_1 MI2_2 MII1_1 MII1_2 MII2_1 MII2_2 PN1_1 PN1_2 PN2_1 PN2_2 CL1_1 CL1_2 CL2_1 CL2_2 MO1_1 MO1_2 MO2_1 MO2_2 IC1_1 IC1_2 IC2_1 IC2_2 TP1_1 TP1_2 TP2_1 TP2_2

# extract 1000000 reads and rename the file
i=1;for file in *.fastq;do name=$(cat my_name|awk '{printf $"'$i'"}'); cat $file|head -4000000 > ../test_data/${name}.fastq; i=`expr $i + 1`; done

# fastQC analysise
fastqc -t 2 -o ../FastQC_result/1_fastqc/ -q *.fastq &

# Remove nextera adapter sequences using trim_galore
trim_galore -q 20 --nextera --paired *.fastq  -o ../trim_report/ &

# rename the file
for file in *.fq;do name=$(ls $file|cut -d_ -f1,2);mv $file trim_${name}.fq;done

# Match up paired-end reads using seqkit pair
for file in *.fq;do name=$(ls $file|cut -d_ -f2,3);cat ${file} | rRNAdust -t 8 ../../U13369.1.fa > ../rRNAdust_result/rRdust_${name};done

# Run mapping through qsub using STAR
qsub qsub.STAR.sh
qsub -I
-------------qsub.STAR.sh-------------
#!/bin/bash
#PBS -q MEDIUM
#PBS -l select=1:ncpus=8
#PBS -V
cd ~/DATA/RNA_seq_1/test_data/seqkit_result/
fa="/user2/tanpaku/mizuguchi/HuYajie/DATA/RNA_seq_1/genomeDIr/hg38_index_dir/"
cpu=8
for file in *.fq;do
        name=$(ls $file|cut -d_ -f1);
        STAR --readFilesIn ${name}_1.paired.fq ${name}_2.paired.fq --runThreadN $cpu --genomeDir $fa --outSAMunmapped Within --outSAMtype BAM SortedByCoordinate --outFileNamePrefix ../STAR_result/${name};
done;

# Count mapped reads for each gene using featureCounts
featureCounts -T 6 -p -t exon -g gene_id -a ~/DATA/RNA_seq_1/genomeDIr/gencode.v38.annotation.gtf -o Counts.txt *.bam