平台：Linux 服务器、Rstudio
数据：NCBI数据库（SRR2932830）
1.下载序列
参考 菜鸟自学之——SRA Toolkit 的下载和使用https://blog.csdn.net/guguaihezi/article/details/81240916
可以从https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/找不同的版本，复制链接，使用wget 链接 在服务器上下载。
wget https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.9.2/sratoolkit.2.9.2-centos_linux64.tar.gz
tar xzvf sratoolkit.2.9.2-centos_linux64.tar.gz
mv sratoolkit.2.9.2-centos_linux64 ~/local/app/ #移动到指定文件夹
cd ~/local/app/ #进入本地程序安装路径
mv sratoolkit.2.9.2-centos_linux64 sratoolkit #去掉版本号是为了避免因升级而需要修改
配置文件（方法一）
vi ~/.bashrc #用vi/vim编辑器修改bashrc文件
i #由command line进入insertion line
export PATH=$PATH:/home/urname/local/app/sratoolkit/bin #改成自己的地址
ESC, :wq #退出vi编辑器并保存文件
source ~/.bashrc #让配置生效
配置文件（方法二）
echo 'export PATH=/home/urname/local/app/sratoolkit/bin:$PATH' >> ~/.bashrc #改成自己的地址
source ~/.bashrc
从NCBI的SRA库里下载数据
prefetch SRR2932830
##或者
nohup prefetch SRR2932830 & #不中断
2.数据质控
（1）将sra文件转化为fastq文件
下载的SRR2932830是一个PE测序结果，所以，需要 --split-files 参数将其分解为两个fastq文件。
如果不加该参数，则只有1个fastq文件（包含了两端测序的结果）
fastq-dump --split-3 SRR2932830 #将双端测序文件拆分为两个fastq文件
##或者
fastq-dump --gzip --split-files SRR2932831 #输出gz的压缩格式
（2）安装fastp或fastqc
安装与配置miniconda3时主要参考 安装fastqc_转录组启动：Miniconda安装|Aspera高速下载测序数据|Fastp质控过滤https://blog.csdn.net/weixin_30166691/article/details/112733596
conda install fastp
conda install fastqc
（3）质控&过滤
参考 RNA-Seq：从fastq到表达矩阵
https://www.jianshu.com/p/3352bfd404f3
fastp (single -end, SE)
fastp -I SRR******.fastq -O SRR******_clean.fastq
fastp (paired -end, PE)
fastp -i Sample1-1_R1.fq.gz -o Sample1-1_R1.clean.fq.gz -I Sample1-1_R2.fq.gz -O Sample1-1_R2.clean.fq.gz
##或者
fastp -i Sample1-1_R1.fq.gz -I Sample1-1_R2.fq.gz -o Sample1-1_R1.clean.fq.gz -O Sample1-1_R2.clean.fq.gz
根据质控结果决定是否需要去接头等
3.比对到参考基因组上——hisat2
（1）下载安装Hisat2
参考 生物软件（一）：生物软件（Hisat2）的安装与运行https://www.jianshu.com/p/5caba78a55a4
进入网站https://daehwankimlab.github.io/hisat2/download/，复制地址
image.png
wget https://cloud.biohpc.swmed.edu/index.php/s/oTtGWbWjaxsQ2Ho/download
unzip hisat2-2.2.1-Linux_x86_64.zip
cd hisat2-2.2.1/
./hisat2
echo 'export PATH~/hisat2-2.2.1:$PATH' >> ~/.bashrc
source ~/.bashrc
（2）下载参考基因组和注释文件
使用了师兄自己组装的参考基因组，其他可参考 RNA-seq的处理流程https://zhuanlan.zhihu.com/p/261360251
（3）建立索引和对比到参考基因组
hisat2-build /public/home/st11/sequence/ReferenceData/hg38/hg38.fa genome #改成自己的.fa文件所在地址
hisat2 -x ./Reference_GRCH38/grch38_index/genome \
-1 ./trim_galore/3b-mc1_R1_val_1.fq.gz \
-2 ./trim_galore/3b-mc1_R2_val_2.fq.gz \
-S 3b-mc1.sam #这里的genome是上一步生成的genome文件位置
4.htseq工具计算count值
（1）下载安装samtools
参考 linux下samtools安装https://blog.csdn.net/Gentlezzx/article/details/121626879
##安装samtools
wget https://github.com/samtools/samtools/releases/download/1.9/samtools-1.9.tar.bz2
tar -xzvf samtools-1.9.tar.bz2
cd samtools-1.9
./configure --prefix=/home/用户名/.local
make
make install
##安装 HTSeq
conda install HTSeq
（2）samtools处理数据
samtools view -@8 -b 3b-mc1.sam > 3b-mc1.bam # sam 转为 bam
samtools sort -n 3b-mc1.bam > 3b-mc1_sort.bam #按照name排序
（3）htseq工具计算count值
htseq-count -f bam -r name -a 10 -t exon -i gene_id -m union\
3b-mc1_sort.bam hg38_HPV16.gtf >30-mc1_counts.txt
（4）ENSG编号转为Symbol
参考 TCGA数据 ENSG编号转为Symbol（基因名称）https://blog.csdn.net/weixin_38987362/article/details/102900370
##打开30-mc1_counts.txt，加上列名，另存为genecounts.csv，在Rstudio中打开
##安装读取相关包
library(stats4)
library(BiocGenerics)
library(parallel)
library("AnnotationDbi")
library("org.Hs.eg.db")
c1g <- read.table("genecounts.csv",header = T,sep = ",",row.names = 1)
gene_up=rownames(c1g)
c1g$symbol <- mapIds(org.Hs.eg.db,keys=gene_up,column="SYMBOL",keytype="ENSEMBL",multiVals="first")
mdata<-rbind(c1g,c1g$symbol)
write.table(mdata,'./mgene.csv',sep=",")
完结撒花~
网友评论