操作记录-2020-09-27：xiaxianyou_RNA_seq

最新推荐文章于 2024-04-17 21:17:24 发布
最新推荐文章于 2024-04-17 21:17:24 发布
阅读量1k 收藏 5
点赞数 1
分类专栏： 操作记录
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/xiaomotong123/article/details/108822735
版权

1. 检查上传数据的完整性

操作记录如下：

#显示各数据的完整性代码
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata$ cat *.txt
ea24c0a81c3b146f52180d03fb387f54  s4301_FKDL202588235-1a_1.clean.fq.gz
44e7699dc34c62f7947e4288ea1f0f72  s4301_FKDL202588235-1a_2.clean.fq.gz
8d30a8af5b877eb4b19803059aab89b7  s43023_FKDL202588236-1a_1.clean.fq.gz
f61f25e52a497519dc864d297a42f3d7  s43023_FKDL202588236-1a_2.clean.fq.gz
1498e251add4acdae5d0055bcb50235e  s43113_FKDL202588238-1a_1.clean.fq.gz
230d90a9cef859cff974a8ca5c160e08  s43113_FKDL202588238-1a_2.clean.fq.gz
250908f71a97358188cf6d9c3c5746fa  s4312_FKDL202588237-1a_1.clean.fq.gz
15551245e7bdeee766d45508712a24df  s4312_FKDL202588237-1a_2.clean.fq.gz
0e2d887e46948f58a7bde596bbf127ec  s43159_FKDL202588239-1a_1.clean.fq.gz
caec36ee9d4cc4022014d69420c1b1a0  s43159_FKDL202588239-1a_2.clean.fq.gz
5b33fdc62ccc9d75ad41dd18145d3215  s43168_FKDL202588240-1a_1.clean.fq.gz
3d6082d763a37856e59ff91c01948349  s43168_FKDL202588240-1a_2.clean.fq.gz
#将各数据的完整代码统一写入待检查的文件中
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata$ echo "ea24c0a81c3b146f52180d03fb387f54  s4301_FKDL202588235-1a_1.clean.fq.gz
> 44e7699dc34c62f7947e4288ea1f0f72  s4301_FKDL202588235-1a_2.clean.fq.gz
> 8d30a8af5b877eb4b19803059aab89b7  s43023_FKDL202588236-1a_1.clean.fq.gz
> f61f25e52a497519dc864d297a42f3d7  s43023_FKDL202588236-1a_2.clean.fq.gz
> 1498e251add4acdae5d0055bcb50235e  s43113_FKDL202588238-1a_1.clean.fq.gz
> 230d90a9cef859cff974a8ca5c160e08  s43113_FKDL202588238-1a_2.clean.fq.gz
> 250908f71a97358188cf6d9c3c5746fa  s4312_FKDL202588237-1a_1.clean.fq.gz
> 15551245e7bdeee766d45508712a24df  s4312_FKDL202588237-1a_2.clean.fq.gz
> 0e2d887e46948f58a7bde596bbf127ec  s43159_FKDL202588239-1a_1.clean.fq.gz
> caec36ee9d4cc4022014d69420c1b1a0  s43159_FKDL202588239-1a_2.clean.fq.gz
> 5b33fdc62ccc9d75ad41dd18145d3215  s43168_FKDL202588240-1a_1.clean.fq.gz
> 3d6082d763a37856e59ff91c01948349  s43168_FKDL202588240-1a_2.clean.fq.gz" >check_md5sum.txt
#执行md5sum 命令对其完整性进行检测
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata$ md5sum -c check_md5sum.txt
s4301_FKDL202588235-1a_1.clean.fq.gz: OK
s4301_FKDL202588235-1a_2.clean.fq.gz: OK
s43023_FKDL202588236-1a_1.clean.fq.gz: OK
s43023_FKDL202588236-1a_2.clean.fq.gz: OK
s43113_FKDL202588238-1a_1.clean.fq.gz: OK
s43113_FKDL202588238-1a_2.clean.fq.gz: OK
s4312_FKDL202588237-1a_1.clean.fq.gz: OK
s4312_FKDL202588237-1a_2.clean.fq.gz: OK
s43159_FKDL202588239-1a_1.clean.fq.gz: OK
s43159_FKDL202588239-1a_2.clean.fq.gz: OK
s43168_FKDL202588240-1a_1.clean.fq.gz: OK
s43168_FKDL202588240-1a_2.clean.fq.gz: OK

2.精简文件名称

操作记录如下：

(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata$ mv s4301_FKDL202588235-1a_1.clean.fq.gz s4301_1.clean.fq.gz
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata$ mv s4301_FKDL202588235-1a_2.clean.fq.gz s4301_2.clean.fq.gz
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata$ mv s43023_FKDL202588236-1a_1.clean.fq.gz s43023_1.clean.fq.gz
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata$ mv s43023_FKDL202588236-1a_2.clean.fq.gz s43023_2.clean.fq.gz
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata$ mv s43113_FKDL202588238-1a_1.clean.fq.gz s43113_1.clean.fq.gz
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata$ mv s43113_FKDL202588238-1a_2.clean.fq.gz s43113_2.clean.fq.gz
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata$ mv s4312_FKDL202588237-1a_1.clean.fq.gz s4312_1.clean.fq.gz
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata$ mv s4312_FKDL202588237-1a_2.clean.fq.gz s4312_2.clean.fq.gz
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata$ mv s43159_FKDL202588239-1a_1.clean.fq.gz s43159_1.clean.fq.gz
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata$ mv s43159_FKDL202588239-1a_2.clean.fq.gz s43159_2.clean.fq.gz
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata$ mv s43168_FKDL202588240-1a_1.clean.fq.gz s43168_1.clean.fq.gz
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata$ mv s43168_FKDL202588240-1a_2.clean.fq.gz s43168_2.clean.fq.gz
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata$ ll
total 16G
drwxrwxr-x 2 zexing zexing  4.0K 9月  27 10:35 .
drwxrwxr-x 4 zexing zexing  4.0K 9月  27 10:12 ..
-rw-rw-r-- 1 zexing zexing   860 9月  27 10:25 check_md5sum.txt
-rw-rw-r-- 1 zexing zexing   142 9月  27 10:09 MD5_s4301_FKDL202588235-1a.txt
-rw-rw-r-- 1 zexing zexing   144 9月  27 10:03 MD5_s43023_FKDL202588236-1a.txt
-rw-rw-r-- 1 zexing zexing   144 9月  27 09:57 MD5_s43113_FKDL202588238-1a.txt
-rw-rw-r-- 1 zexing zexing   142 9月  27 09:51 MD5_s4312_FKDL202588237-1a.txt
-rw-rw-r-- 1 zexing zexing   144 9月  27 09:46 MD5_s43159_FKDL202588239-1a.txt
-rw-rw-r-- 1 zexing zexing   144 9月  27 09:39 MD5_s43168_FKDL202588240-1a.txt
-rw-rw-r-- 1 zexing zexing  1.3G 9月  27 10:09 s4301_1.clean.fq.gz
-rw-rw-r-- 1 zexing zexing  1.3G 9月  27 10:06 s4301_2.clean.fq.gz
-rw-rw-r-- 1 zexing zexing  1.5G 9月  27 10:03 s43023_1.clean.fq.gz
-rw-rw-r-- 1 zexing zexing  1.5G 9月  27 10:00 s43023_2.clean.fq.gz
-rw-rw-r-- 1 zexing zexing  1.5G 9月  27 09:57 s43113_1.clean.fq.gz
-rw-rw-r-- 1 zexing zexing  1.6G 9月  27 09:54 s43113_2.clean.fq.gz
-rw-rw-r-- 1 zexing zexing  1.2G 9月  27 09:51 s4312_1.clean.fq.gz
-rw-rw-r-- 1 zexing zexing  1.2G 9月  27 09:48 s4312_2.clean.fq.gz
-rw-rw-r-- 1 zexing zexing  1.5G 9月  27 09:46 s43159_1.clean.fq.gz
-rw-rw-r-- 1 zexing zexing  1.6G 9月  27 09:42 s43159_2.clean.fq.gz
-rw-rw-r-- 1 zexing zexing  948M 9月  27 09:39 s43168_1.clean.fq.gz
-rw-rw-r-- 1 zexing zexing 1017M 9月  27 09:37 s43168_2.clean.fq.gz

3. 使用FastQC软件对数据进行质控

操作记录如下：

#编辑脚本如下：
#上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program:
#       This program is used for analysis of RNA-seq data by FastQC.
#History:
# 2020/09/27         zexing              First release
#fastqc命令为质控命令
#Usage: fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN
#简写代码：fastqc -t 8 -o <out-dir> seqfile1
#调用程序fastqc，参数-t设置线程数为8，参数-o设置结果输出的目录，参数-c可以加入污染物选项（接头信息），最后为读入文件。
fastqc -t 16 -o /f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/fastqc_report/ /f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/*.fq.gz

#后台运行命令如下：
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/scripts_log$ nohup bash fastqc_script &
[1] 123084

4. 使用Hisat2软件对测序结果进行比对

操作记录如下：

#编辑脚本如下：
#! /bin/bash
#上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program:
#       This program is used for aligning of RNA-seq data.
#History:
# 2020/09/27         zexing              First release
#对变量${i}利用 for ${i} in A B C D 的方式遍历指定
for i in s4301 s4312 s43023 s43113 s43159 s43168
do
#hisat2命令为序列比对命令
#Usage: hisat2 [options]* -x <ht2-idx> {-1 <m1> -2 <m2> | -U <r> | --sra-acc <SRA accession number>} [-S <sam>]
#简写代码：hisat2 -t -p 8 -x <ht2-idx> -1 fq.gz -2 fq.gz -S <sam>
#调用程序hisat2，参数-t显示时间，参数-p设定线程数，参数-x指定参考基因组索引文件的前缀（目录及文件前缀）
#参数-1指定双向测序的第一条.fq.gz测序数据，如果多组数据，使用逗号将文件分隔
#参数-2指定双向测序的第二条.fq.gz测序数据，如果多组数据，使用逗号将文件分隔
#参数-S指定比对结果的输出目录及名称
hisat2 -t -p 16 -x /f/xudonglab/zexing/reference/UCSC_mm10/hisat2_index/hisat2_index_mm10 \
-1 /f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/${i}_1.clean.fq.gz \
-2 /f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/${i}_2.clean.fq.gz \
-S /f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/aligned/sam/${i}.sam
done

#后台运行命令如下：
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/scripts_log$ nohup bash Hisat2_script &
[1] 123564

5. 使用SAMtools软件对文件进行格式转换、排序

操作记录如下：

#编辑脚本如下：
#!/bin/bash
#上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program:
#这里写一个脚本，将samtools的view、sort、index和flagstat四个命令串联在一起，对RNA-seq数据批量处理。
#History:
# 2020/09/27         李泽兴      First release
#对变量${i}利用 for ${i} in A B C D 的方式遍历指定
#samtools view为格式转换命令
#Usage: samtools view [options] <in.bam> | <in.sam> | <in.cram> [region...]
#-@设置线程数为8， -S输入sam文件，-1b使用快速压缩比生成bam文件，-o设定标准输出文件名。
#简写代码：samtools view -@ 8 -S {$i}.sam -1b -o {$i}.bam
#samtools sort为排序命令
#Usage: samtools sort [options..] [in.bam]
#-@设置线程数为8，-l设置压缩比为5，-o设定标准输出文件名，最后一个为输入.bam文件，该命令默认为染色体位置排序，不影响后期操作。
#简写代码：samtools sort -@ 8 -l 5 -o {$i}.bam.sort {$i}.bam
#index为建立索引命令，目前还不知道这样建立索引有何用处，暂且写下来以备后续使用。
#Usage: samtools index [-bc] [-m INT] <in.bam> <out.index>
#-@设置线程数为8，<in.bam>输入bam文件，<out.index>输出index文件。
#简写代码：samtools index -@ 8 {$i}.bam.sort {$i}.bam.index
#flagstat为查看比对情况命令，目前暂且对这部分不是很熟悉，暂且写下来以备后续使用。
#Usage: samtools flagstat [options] <in.bam>
#-@设置线程数为8，对于输出内容重定向写入>指定文件
#简写代码：samtools flagstat -@ 8 {$i}.bam.sort > {$i}.bam.flagstat
for i in s4301 s4312 s43023 s43113 s43159 s43168
do
samtools view -@ 16 -S /f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/aligned/sam/${i}.sam \
-1b -o /f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/aligned/bam/${i}.bam
samtools sort -@ 16 -l 5 -o /f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/aligned/bam.sort/${i}.bam.sort \
/f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/aligned/bam/${i}.bam
samtools index -@ 16 /f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/aligned/bam.sort/${i}.bam.sort \
/f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/aligned/bam.index/${i}.bam.index
samtools flagstat -@ 16 /f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/aligned/bam.sort/${i}.bam.sort \
> /f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/scripts_log/${i}.bam.flagstat
done

#后台运行命令如下：
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/scripts_log$ nohup bash SAMtools_script &
[1] 127169

6. 使用RSeQC软件对比对结果进行质控

操作记录如下：

#编辑脚本如下：
#!/bin/bash
#上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
#program:
#  This program is running RSeQC software command and saving the output results to special files.
#History:
#    2020/09/27            zexing         First release
#对变量${i}利用 for ${i} in A B C D 的方式遍历指定
#bam_stat.py命令为python一个程序，用于查看比对总体情况。
#Usage: bam_stat.py [options]
#-i参数设定输入文件，重定向写入>指定结果输出文件。
#简写代码：bam_stat.py -i ${i}.bam.sort > ${i}_bam_stat.log
#read_distribution.py命令为python一个程序，用于查看基因组覆盖率。
#Usage: read_distribution.py [options]
#-i参数设定输入文件，-r参数指定参考的bed文件，重定向写入>指定结果输出文件。
#简写代码：read_distribution.py -i ${i}.bam.sort -r RefSeq.bed > ${i}_distribution.log
for i in s4301 s4312 s43023 s43113 s43159 s43168
do
bam_stat.py -i /f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/aligned/bam.sort/${i}.bam.sort \
> /f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/scripts_log/${i}_bam_stat.log
read_distribution.py -i /f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/aligned/bam.sort/${i}.bam.sort \
-r /f/xudonglab/zexing/reference/BED_file/mm10_RefSeq.bed \
> /f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/scripts_log/${i}_distribution.log
done

#后台运行命令如下：
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/scripts_log$ nohup bash QSeQC_script &
[1] 147889

7. 使用Stringtie软件进行拼接和定量

操作记录如下：

#编辑脚本如下：
#!/bin/bash
#上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
#program
#     This program is to perform HTSeq-count assay for RNA-seq data.
#History
#     2020/09/27       zexing            First release
#对变量${i}利用 for ${i} in A B C D 的方式遍历指定
#创建针对不同样品的ballgown/"$i"目录，以存放输出结果
#stringtie为转录本拼接和定量的软件
#利用stringtie分别输出gene count数目和FPKM值，gene count数用户基因差异分析，FPKM值为基因表达值，用于heatmap等绘制
#Usage: stringtie <aligned_reads.bam> [options]
#aligned_reads.bam 是输入文件，该输入文件要求必须按基因组位置排序
#参数-o [<path/>]<out.gtf>，设置StringTie组装转录本的输出GTF文件的路径和文件名。
#参数-p <int>指定组装转录本的线程数（CPU）
#参数-G <ref_ann.gff>使用参考注释基因文件指导组装过程，格式GTF/GFF3。输出文件中既包含已知表达的转录本，也包含新的转录本。
#参数-B 应用该选项，则会输出Ballgown输入表文件（* .ctab），其中包含用-G选项给出的参考转录本的覆盖率数据。
#参数-e 限制reads比对的处理，仅估计和输出与用-G选项给出的参考转录本匹配的组装转录本。使用该选项，则会跳过处理与参考转录本不匹配的组装转录本，这将大大的提升了处理速度。
#参数-A <gene_abund.tab> 输出基因丰度的文件（制表符分隔格式）
#简写代码：stringtie {$i}.bam.sort -o path/name.gtf -p 16 -G gene.gtf -eB -A path/name.tab
#此次主要输出两个文件：gtf文件用于转换gene count文件；tab文件包含FPKM数值。
for i in s4301 s4312 s43023 s43113 s43159 s43168
do
mkdir /f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/aligned/ballgown/"$i"
stringtie /f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/aligned/bam.sort/"$i".bam.sort \
-o /f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/aligned/ballgown/"$i"/"$i".gtf \
-p 16 -G /f/xudonglab/zexing/reference/UCSC_mm10/mm10_genes.gtf -e -B \
-A /f/xudonglab/zexing/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/aligned/ballgown/"$i"/"$i".gene.tab
done

#后台运行命令如下：
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/scripts_log$ nohup bash Stringtie_script &
[2] 148677

8. 使用prepDE.py脚本提取read counts数值

注意：prepDE.py为python2的脚本，应该先将python设置为低版本后再运行脚本。

操作记录如下：

#从之前的文件夹中将prepDE.py脚本拷贝至当前文件夹中的ballgown子文件夹中
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/aligned/ballgown$ cp /f/xudonglab/zexing/projects/zhaoxiujuan/aligned/ballgown_1/prepDE.py ./

#退出当前conda环境
(base) zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/aligned$ conda deactivate

#检查服务器中的python版本信息
zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/aligned$ python --version
Python 2.7.12

#使用python命令直接运行脚本
zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/aligned/ballgown$ python prepDE.py

#查看运行结果
zexing@DNA:~/projects/xiaxianyou/RNA_seq/2020_09_23/cleandata/aligned/ballgown$ ll
total 2.0M
drwxrwxr-x 8 zexing zexing 4.0K 9月  28 17:02 .
drwxrwxr-x 8 zexing zexing 4.0K 9月  28 17:01 ..
-rw-rw-r-- 1 zexing zexing 860K 9月  28 17:02 gene_count_matrix.csv
-rw-rw-r-- 1 zexing zexing  12K 9月  28 16:59 prepDE.py
drwxrwxr-x 2 zexing zexing 4.0K 9月  27 15:25 s4301
drwxrwxr-x 2 zexing zexing 4.0K 9月  27 15:12 s43023
drwxrwxr-x 2 zexing zexing 4.0K 9月  27 15:14 s43113
drwxrwxr-x 2 zexing zexing 4.0K 9月  27 15:09 s4312
drwxrwxr-x 2 zexing zexing 4.0K 9月  27 15:17 s43159
drwxrwxr-x 2 zexing zexing 4.0K 9月  27 15:18 s43168
-rw-rw-r-- 1 zexing zexing 1.1M 9月  28 17:02 transcript_count_matrix9以下.csv
#其中"gene_count_matrix.csv"即是DESeq2的输入文件。

以下分析在本地机的Rstudio中完成

9.利用DESeq2包进行差异基因分析

代码如下：

#This script is used for analysis of xiaxianyou RNA-seq data
#History
# Lizexing           2020-09-28             First release
#genecount文件来源于Stringtie软件分析，后面为本地电脑操作
# DESeq2进行差异分析 ----------------------------------------------------------------

#清空环境变量
rm(list=ls())
#设置工作目录
setwd("G:/xiaxianyou/RNA-seq/2020_09_28/gene_count/")
#读入基因表达值，设定行名为gene_id
gene_count <- read.csv("gene_count_matrix.csv",stringsAsFactors = F)
#对gene_id一列进行拆分，去除重复名称
library(stringr)
#设置空的"gene_count_1"向量，行数与上面的测序结果一致
gene_count_1<-rep(NA,nrow(gene_count))
#利用for循环，对gene_count数据框中的重复列进行拆分提取
for (i in 1:nrow(gene_count)){
   
  gene_count_1[i] <- unlist(str_split(gene_count[i,1], pattern = "\\|"))[1]
}
#显示拆分后的结果
head(gene_count_1)
#对原数据框中的特定序列重新赋值
gene_count$gene_id <- gene_count_1
#显示文件的前6行信息
head(gene_count)
#将第一列作为文件的行名
rownames(gene_count) <- gene_count[,1]
gene_count <-gene_count[,-1]
#显示文件的前6行信息
head(gene_count)
#将各组数据分开
gene_count_group_1 <- gene_count[, 1:4]
gene_count_group_2 <- gene_count[, c(1, 2, 5, 6)]
#将该文件保存至对应目录
write.csv(gene_count_group_1, file = "G:/xiaxianyou/RNA-seq/2020_09_28/gene_count/gene_count_group_1.csv", row.names = TRUE)
write.csv(gene_count_group_2, file = "G:/xiaxianyou/RNA-seq/2020_09_28/gene_count/gene_count_group_2.csv", row.names = TRUE)

#加载DESeq2包
library(DESeq2)
#DESeq2需要三种矩阵，分别为countData表达矩阵,colData样品信息矩阵及design差异表达矩阵
#countData为表达矩阵即gene_count
#colData为样品信息矩阵即coldata
#design为差异表达矩阵即批次和条件（对照、处理）等
#设置condition样品组别、重复数
condition_group_1 <- factor(c(rep("scr", 2), rep("shA17_B", 2)), levels = c("scr","shA17_B"))
condition_group_2 <- factor(c(rep("scr", 2), rep("shA17_C", 2)), levels = c("scr","shA17_C"))
#显示condition设置信息
condition_group_1
condition_group_2
#设置group组对应的样品信息矩阵colData
colData_group_1 <- data.frame(row.names = colnames(gene_count_group_1), condition_group_1)
colData_group_2 <- data.frame(row.names = colnames(gene_count_group_2), condition_group_2)
#显示colData设置信息
colData_group_1
colData_group_2
#在R里面用于构建公式对象，~左边为因变量，右边为自变量。
#标准流程：dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design= ~ batch + condition) 
#countData为表达矩阵即countdata
#colData为样品信息矩阵即coldata
#design为差异表达矩阵即批次和条件（对照、处理）等
#对dds_group进行标准流程构建
dds_group_1 <- DESeqDataSetFromMatrix(gene_count_group_1, colData_group_1, design = ~condition_group_1)
dds_group_2 <- DESeqDataSetFromMatrix(gene_count_group_2, colData_group_2, design = ~condition_group_2)
#对原始dds_group进行normalize
dds_group_1 <- DESeq(dds_group_1)
dds_group_2 <- DESeq(dds_group_2)
#显示dds信息
dds_group_1
dds_group_2

# 对差异分析结果进行保存 -------------------------------------------------------------

#使用DESeq2包中的results()函数，提取差异分析的结果
#Usage:results(object, contrast, name, .....）
#将提取的差异分析结果定义为变量"res" 
#contrast: 定义谁和谁比较，处理组在前，对照组在后
#将group组提取分析结果并保存为res
res_group_1 = results(dds_group_1, contrast=c("condition_group_1","shA17_B","scr"))
res_group_2 = results(dds_group_2, contrast=c("condition_group_2","shA17_C","scr"))

#对结果res利用order()函数按pvalue值进行排序
#创建矩阵时，X[i,]指矩阵X中的第i行，X[,j]指矩阵X中的第j列
#order()函数先对数值排序，然后返回排序后各数值的索引，常用用法：V[order(V)]或者df[order(df$variable),]
#对res_group组进行排序
res_group_1 = res_group_1[order(res_group_1$pvalue),]
res_group_2 = res_group_2[order(res_group_2$pvalue),]

#显示res结果前6行信息
head(res_group_1)
head(res_group_2)

#对res_group矩阵进行总结，利用summary命令统计显示一共多少个genes上调和下调
summary(res_group_1
最低0.47元/天 解锁文章
垚垚爸爱学习
关注
  • 1
    点赞
  • 踩
  • 5
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
Video-Cap::clapper_board:视频字幕
05-16
视频字幕 此存储库是ICCV '15论文-Tensorflow 1.0中的的实现 表现 方法 BLEU @ 1分数 seq2seq * 0.28 * seq2seq是纸张模型的复制品 ... 您可以使用“ data / test_features”目录来简化操作 训练密码 pip insta
2022-TCGA数据库重大更新后RNASeq的STAR-Counts数据的下载与整理
BioInfoNotes
04-06 1万+
TCGA|GEO|文献阅读|数据库|理论知识R语言|Bioconductor| 服务器与Linux最近有粉丝留言，TCGA数据库发生更新，下载的数据和之前的不一样。比如转录组，之前是HTSeq流程的数据，现在是STAR-Counts的数据。具体的数据信息参考：https://docs.gdc.cancer.gov/Data/Release_Notes...
利用BLAST进行序列比对和寻找同源基因
热门推荐
dpsh_6的博客
03-10 3万+
NCBI BLAST是最常用的短序列局部比对软件之一，学习一下软件参数，并这记录一下我的软件操作的命令行。 NR 库分割参考Gu Kai老师的笔记：http://www.bioinfo-scrounger.com/archives/673 1. 软件版本 我部署在服务器上的是version 2.8.1 ，因为其支持物种注释以及分割的NR子库的索引作为比对的库，使用比较方便，但...
基因驱动 CovEx 代码阅读（一）
GOD_胖慧的博客
07-06 957
CovEx 源代码　https://sourceforge.net/projects/cancer-pathway/ 首先重点阅读 README1.txt 一步一步执行 安装 Gurobi Optimizer 注意： 这里需要校园网环境 官网上下载到gurobi8.0.1_linux64.tar.gz 安装产生gurobi.lic 具体步骤参考： https://blog.csdn.net/s...
AC自动机+DP HDOJ 2457 DNA repair（DNA修复）
weixin_33736649的博客
07-14 142
题目链接 题意： 　　给n串有疾病的DNA序列，现有一串DNA序列，问最少修改几个DNA，能使新的DNA序列不含有疾病的DNA序列。 思路： 　　构建AC自动机，设定end结点，dp[i][j]表示长度i的前缀串走到自动机的j点最少需要修改几个DNA。状态转移方程。那么只要转移到下一个的DNA不是end结点就能转移，如果下一个DNA不和原序列不一样就+1。 #include &...
Stringtie 计算转录组的 Raw Counts
公众号/简说基因，知乎/简宝玉
03-20 1944
Stringtie 自带一个脚本prepDE.py用于计算转录组的 Raw Counts，用法如下：Usage: prepDE.py [options] Generates two CS...
rnaseq:RNA-seq分析
04-28
RNA-seq管线 设置 在raw运行中，先执行./fetch_data.sh ，然后执行./create_index.sh 。 在src运行./fetch_binaries.sh （需要的MacOS或Linux）， ./setup_samtools.sh和./setup_rsem.sh （需要各种库编译）。 ...
grape-nf:使用Nextflow的自动RNA-seq管道
05-19
Grape提供了广泛的RNA-Seq分析渠道。 它允许创建自动化和集成的工作流程来管理和分析RNA-Seq数据。 它使用作为执行后端。 请查看以获取更多信息。 在IHEC联盟中，Grape已被用于RNA-seq整合分析。 根据IHEC建议，...
kallisto:接近最佳的RNA-Seq定量
04-28
在具有标准RNA-Seq数据的基准测试中， kallisto可以在Mac台式计算机上不到3分钟的时间内，仅使用读取序列和本身需要10分钟即可完成的转录组索引，就可以定量分析3000万个人的大量RNA-seq读数。 读数的伪比对保留了...
SARTools：RNA-Seq工具的统计分析
02-05
RNA-Seq技术是一种广泛应用的高通量测序方法，用于研究细胞或组织中的转录组。SARTools（Statistical Analysis for RNA-Seq Tools）是一个专门针对RNA-Seq数据进行统计分析的工具集，旨在帮助生物信息学家和研究人员...
操作记录-2020-11-08：精简代码处理RNA_seq数据
垚垚爸的博客
11-08 2831
今天准备尝试编写一组精简代码用于处理RNA_seq数据，希望能成功吧。 1.建立相应目录 对新数据建立对应实验人员（lizexing）、测序类型（RNA_seq）和日期（2020_10_20）的目录。 # 建立后如下： (base) zexing@DNA:~/projects/lizexing/RNA_seq/2020_11_08$ # 新建对应的各级目录 mkdir raw_data clean_data ballgown bam bam_sort sam fastqc_report GSEA MD5_
webbrowser1 脚本报错_HISAT2+Stringtie+DESeq2 （prepDE.py脚本报错向)
weixin_29817925的博客
12-20 1011
2019/10/14 原作者出了解决方案，更新了prepDE.py脚本。=====================原文==============================写这篇文章的主要原因还是出于这两天来找一个很小的bug找不到，最后发现版本不兼容的故事。背景：HISAT2+Stringtie+Ballgown本来是一组黄金组合，但是由于我的生物学重复（biological replic...
转录组上游分析，Count计算
最新发布
BioinfoDu
04-17 1414
featureCounts是subread中脚本，可以使用subread流程进行定量，在这里直接使用前面mapped的bam文件进行转录本定量。，主要发表或收录生物信息学的教程，以及基于R的分析和可视化（包括数据分析，图形绘制等）；因此，个人建议仍是获得count值，后再进一步的分析，这样的方法更利于下游分析。是stringtie软件中自带的获得转录本丰度的脚本，6.1中的。文件进行转录本表达量的比对，获得转录本的FPKM，此后使用。在网络中吗，都有比较详细的教程，大家可以自己去学习。
操作记录-2022-09-12: xieruyu_project
垚垚爸的博客
09-12 707
鸡的转录基因组分析
2021-07-14:sunruiqi_RNA_seq
垚垚爸的博客
07-15 503
今天准备尝试编写一组精简代码用于处理RNA_seq数据，希望能成功吧。 1.建立相应目录 对新数据建立对应实验人员（sunruiqi）、测序类型（RNA_seq）和日期（2021_07_14）的目录。 # 建立后如下： (base) zexing@bio:~/projects/sunruiqi/RNA_seq/2021_07_14$ # 新建对应的目录 mkdir raw_data clean_data ballgown bam bam_sort sam fastqc_report GSEA MD5_tx
R语言——GO分析和KEGG分析
codyliu_bo的博客
10-30 1720
进行GO和KEGG分析
生物信息学——常见的OrgDB包
qq_45794091的博客
06-20 1812
物种 OrgDB 人（Humanm) org.Hs.eg.db 小鼠（Mouse） org.Mm.eg.db 猪（Pig) org.Ss.eg.db 按蚊（Anopheles） org.Ag.eg.db 拟南芥（Arabidopsis） org.At.tair.db 牛（Brovine) org.Bt.eg.db 蠕虫（Worm） org.Ce.eg.db 犬（Canine） org.Cf.eg.db 苍蝇（Fly) org.Dm.eg.db 斑马鱼（Ze
Linux的RNA-seq分析（stringtie-FPKM和HTseq计数）
weixin_43927366的博客
03-15 1569
确保安装了stringtie如果要研究新的转录本，则需要完成所有样本的转录本的组装并整合在stingtie-merged.gtf中，详细见Linux系统RNA-seq分析（stringtie转录本组装)。新建SRR3418005_FPKM子文件夹，并在其中打开终端：这个选项告诉stringtie不仅要输出组装得到的转录本，还要输出那些表达量较低、可能不被认为是完整转录本的“额外”转录本片段。这有助于捕获那些可能只部分被测序覆盖的转录本区域。-p 2：这个选项指定了stringtie。
2020-09-30
qq_40310710的博客
09-30 151
9.18 新氧深度学习岗位： 编程问题： python： Numpy():np.where(mask>0) Pandas(): Pytorch(): 继承torch.nn 网络：超分辨率 GAN：用GAN做超分辨率：GAN里面的损失函数感知损失、对抗损失 ResNet：主要解决的问题：梯度消失和梯度爆炸、对于深层网络存在的退化问题 BatchNormalization：也是用来解决梯度消失和爆炸问题的，但不能解决模型退化问题。 感受野：两个33的感受野和55的感受野大小相同 前者的计算量较小 准确率、

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值