使用SnapATAC分析单细胞ATAC-seq数据（一）: SnapATAC简介与安装

SnapATAC (Single Nucleus Analysis Pipeline for ATAC-seq) 是一个能够快速、准确和全面分析单细胞ATAC-seq数据的R包，它可以对单细胞ATAC-seq数据进行常规的数据降维、聚类和批次校正分析，鉴定远端调控元件并预测其调控的靶基因，调用chromVAR软件进行motif分析，同时还可以将scRNA-seq和scATAC-seq数据进行整合分析等。

Rongxin Fang, Sebastian Preissl, Xiaomeng Hou, Jacinta Lucero, Xinxin Wang, Amir Motamedi, Andrew K Shiau, Eran A Mukamel, Yanxiao Zhang, M Margarita Behrens, Joseph Ecker, Bing Ren Fast and Accurate Clustering of Single Cell Epigenomes Reveals Cis-Regulatory Elements in Rare Cell Types bioRxiv 615179; doi: >(一)安装bowtie
Bowtie可以在个人计算机上使用，也可以在CSC服务器上使用终端连接。请参阅以下文档的第一部分，了解如何在笔记本电脑上安装Bowtie。特别是对他们的计算机没有管理员权限的那些应该确保软件的正确安装和功能。Bowtie也可以在服务器计算机上远程使用。我们将提供临时帐户访问CSC，但你将需要一个安全Shell终端程序进行通信。默认情况下，Mac和Linux上都有这样的程序，但需要安装Windows。普遍的实现是PuTTY。即使终端程序不用于读取映射，也将需要其他练习，并且应该可用。Bowtie的安装：从下载页面下载相应的版本(Linux，Mac或Win,小编使用的是在Linux下进行)。将zip文件解压缩到新的目录中，并转到该目录。下载的bowtie包装包含大肠杆菌基因组的预先建立的指数，以及从该基因组模拟的一组1000个35bp的读数。要使用Bowtie对齐这些读取，请键入以下命令。bowtiee_colireads/e_coli_1000fqmap_resulttxt
如果你收到错误消息"commandnotfound"，请尝试在"bowtie"(/bowtie)之前添加"/"。
(二)使用Bowtie
(1)Mapping
要使用Bowtie对齐示例读取，请发出以下命令。bowtiee_colireads/e_coli_1000fqmap_resulttxt
如果你收到错误消息"commandnotfound"，请尝试在"bowtie"(/bowtie)之前添加"/"。"e_coli"与"indexes/e_coli"相同。你可以在文本编辑器中打开map_resulttxt。每行都是一个读取对齐。对齐读取的名称显示在第一列中。对于Mac和Linux，使用"少"会更好。
lessmap_resulttxt#extrareading
ReadthemanualinthefolderorwebsitetogetadeeperunderstandinghowBowtieworksandfurtheroptionsinBowtie
我们来看看Bowtie在1中使用的一些不同的选项，报告所有有效的对齐方式与一些不匹配。
/bowtie-a-v2e_coli--suppress1,5,6,7-cATGCATCATGCGCCAT-a/--all报告每个读取或对的所有有效对齐(默认值：off)
-v
最多不相匹配的报告对齐
-c
查询序列在命令行
--suppress
上以默认输出模式抑制输出列
2限制对齐
$/bowtie-k3-v2e_coli--suppress1,5,6,7-cATGCATCATGCGCCAT-k
每次读取或配对时报告有效对齐(默认值：1)。
3不匹配排名
$/bowtie-a--best-v2e_coli--suppress1,5,6,7-cATGCATCATGCGCCAT
所有相同的对齐方式按最佳到最坏的顺序进行报告
4只有最不匹配
$/bowtie-a--best--strata-v2--suppress1,5,6,7e_coli-cATGCATCATG
(2)配对对齐
当使用-1和-2选项指定正确配对的读取文件时，Bowtie可以对齐配对端读取(对于原始，FASTA或FASTQ读取文件)
/bowtiee_coli-1reads/e_coli_1000_1fq-2reads/e_coli_1000_2fqmap_pairedtxt
SAMtools(>===============建立索引===========================

bismark_genome_preparation --bowtie2 / --verbose  index/

================比对=============================

bismark --bowtie2 -N 0 -L 20 --quiet --un --ambiguous --sam --nucleotide_coverage --genome index/ --samtools_path /gpfs03/home/jingjing/software/samtools-19/ -o SRR10401142 -1 /fastq/SRR104011421_1fastqgz -2 /fastq/SRR104011421_2fastqgz

================去重===============================

deduplicate_bismark --samtools_path /gpfs03/home/jingjing/software/samtools-19/ -p SRR104011421_1_bismark_bt2_pesam --output_dir dedup

=====================提取甲基化信息==================

bismark_methylation_extractor -p --comprehensive --no_overlap --bedGraph --counts --buffer_size 200G --report --cytosine_report --samtools_path /gpfs03/home/jingjing/software/samtools-19/ --genome_folder index/ dedup/SRR104011421_1_bismark_bt2_pededuplicatedbam -o SRR10401142/

生成处理报告：

bismark2report

它包括了比对信息，甲基化信息，M-bias等，可以对数据有一个大概的认知。
结果合并正反链的数据后会输出CpG/CHG/CHH三种类型的甲基化文件，包含了胞嘧啶所有的组合形式，但实际上我们自然最关注的是CpG位点的甲基化。其中

CpG_context_SRR104011421_1_bismark_bt2_pededuplicatedtxt即CpG甲基化位点的文件。
# 第一列为测序信息

# 第二列为甲基化状态 + 代表甲基化 -代表未甲基化

# 第三列代表chromosome

# 第四列代表location

# 第五列代表methylation call，简单来说大写的就是甲基化的（因为还有CHG，CHH的数据，分别对应x, X , h, H）

SRR104011421_1_bismark_bt2_pededuplicatedbismarkcovgz文件则给了每个位点的甲基化比例，为下一步确定CpG岛提供了基础，其数据形式如下：

其中：# 第一列代表chromosome  # 第二，三列代表location  # 第四列代表甲基化百分比

# 第五列代表甲基化数目  # 第六列代表未甲基化数目

1 比对的是：使用idba_ud拼接的AER314-4raw_data基因组与转录组数据。

2 bowtie2做index（bowtie2使用conda安装）
建索引：bowtie2-build AER314-4_scaffoldfa AER314-4_scaffoldfa
3 reads mapping到参考基因组——tophat2软件：基于bowtie2(tophat安装见软件安装）

命令：tophat2 -p 12 -o AER314-4_output /home/test04/lyr/rna-seq/02align_out/AER314-4_scaffoldfa /home/test04/lyr/rna-seq/01data/YSH-qurRNA-42-314-4_L001_R1fastq /home/test04/lyr/rna-seq/01data/YSH-qurRNA-42-314-4_L001_R2fastq

4 然后就很顺利的跑出来结果了

使用公司服务器，12个线程，大概五个小时就跑完啦。
5 cufflink

[   Cufflinks输出结果

cufflinks的输入文件是sam或bam格式。并且sam或bam格式的文件必须排好序。（The SAM file supplied to Cufflinks must be sorted by          reference position）Tophat的输出结果sam或bam已经排好了序。针对其他的未排序的sam或bam文件采用如下排序方式：

sort -k 3,3 -k 4,4n hitssam > hitssamsorted

1 transcriptsgtf

该文件包含Cufflinks的组装结果isoforms。前7列为标准的GTF格式，最后一列为attributes。其每一列的意义：

列数  列的名称  例子        描述

1    序列名    chrX        染色体或contig名; 2    来源      Cufflinks  产生该文件的程序名; 3    类型      exon        记录的类型，一般是transcript或exon; 4    起始      1          1-base的值; 5    结束      1000        结束位置; 6    得分      1000        ; 7    链        +          Cufflinks猜测isoform来自参考序列的那一条链，一般是'+','-'或'';8    frame              Cufflinks不去预测起始或终止密码子框的位置; 9    attributes        详见下

每一个GTF记录包含如下attributes：

Attribute      例子      描述

gene_idCUFF1Cufflinks的gene id;transcript_idCUFF11  Cufflinks的转录子 id; FPKM          101267  isoform水平上的丰度, F ragments P er K ilobase of exon model per M illion mapped fragments; frac          07647    保留着的一项，忽略即可，以后可能会取消这个;conf_lo        007      isoform丰度的95%置信区间的下边界，即 下边界值 = FPKM ( 10 - conf_lo );conf_hi        01102    isoform丰度的95%置信区间的上边界，即 上边界值 = FPKM ( 10 + conf_hi ); cov            100765 计算整个transcript上read的覆盖度;full_read_support  yes  当使用 RABT assembly 时，该选项报告所有的introns和exons是否完全被reads所覆盖

2 ispformsfpkm_tracking

isoforms(可以理解为gene的各个外显子)的fpkm计算结果

3 genesfpkm_tracking

gene的fpkm计算结果Cuffmerge简介

Cuffmerge将各个Cufflinks生成的transcriptsgtf文件融合称为一个更加全面的transcripts注释结果文件mergedgtf。以利于用Cuffdiff来分析基因差异表达。

2 使用方法

$ cuffmerge [options]

输入文件为一个文本文件，是包含着GTF文件路径的list。常用例子：

$ cuffmerge -o /merged_asm -p 8 assembly_listtxt

3 使用参数

-h | --help

-o  default: /merged_asm

将结果输出至该文件夹。

-g | --ref-gtf将该reference GTF一起融合到最终结果中。

-p | --num-threads  defautl: 1

使用的CPU线程数

-s | --ref-sequence /该参数指向基因组DNA序列。如果是一个文件夹，则每个contig则是一个fasta文件；如果是一个fasta文件，则所有的contigs都需要在里面。Cuffmerge将使用该ref-sequence来帮助对transfrags分类，并排除repeats。比如transcripts包含一些小写碱基的将归类到repeats  ]

cufflinks：

<1>命令：cufflinks -p 4 -o test_cuff /home/andengdi/lyr/rna-seq/02-align_out/test_output/accepted_hitsbam

流程及结果

5 用相同的方法将其他两个样本跑一遍。

sisichen �
关注
国家基因组科学数据中心（NGDC）---组学原始数据如何上传GSA 原创
2022-04-25 14:44:31
sisichen �
码龄4年
关注
文章目录
前言
一、什么是NGDC？
二、NGDC的发展历程
三、什么是GSA？
四、为什么选择上传数据到GSA？
五、如何上传测序原始数据至GSA？（重点！！附详细步骤！！）
1 准备要上传的数据
2 计算MD5码
3进入NGDC主页，登入账户
4 填写数据信息
第一步：建立Bioproject。
第二步：建立BioSample。
第三步：创建GSA。
进入GSA数据库
新建GSA
填写信息
下载表格文件
5 数据上传：
(1) 通过FTP软件 上传（上传需要流量！！如果小数据可以用）
(2) 通过服务器上传（推荐！！）:如果实验室有服务器的话，推荐服务器上传，步骤如下：（服务器上要先安装ftp ）
（3）邮寄硬盘
6等待审核
总结
前言
在发表文章之前我们需要将测序的原始数据上传到一个公共库，并在文中提供accession number，实现数据的公开共享，这是国际惯例。以前我们上传数据时只能上传到美国国立生物技术信息中心（NCBI）、欧洲生物信息学研究所（EBI）、日本核酸数据库(DDBJ)，现在中国科学院北京基因组研究所（国家生物信息中心）国家基因组科学数据中心 (CNCB-NGDC)—中国的 “NCBI” 已经建立并日渐完善。组学原始数据归档库（GSA）是组学原始数据汇交、存储、管理与共享系统，是国内首个被国际期刊认可的组学数据发布平台。GSA已获得多个国际期刊认可，并已被国际著名出版商Elsevier收录为指定的基因数据归档库，其权威性得到国内外100余家学术杂志的认可。GSA已通过FAIRsharing认证，获得Wiley出版集团认可，因此我们不用担心上传数据到GSA不被期刊认可，也不用再舍近求远上传数据到NCBI，作为中国人，我们一定要支持我们NGDC中的数据库。本文介绍了如何上传测序原始数据到GSA，附详细 *** 作步骤。
一、什么是NGDC？
国家基因组科学数据中心（>

Fastqc
Fastqc website ( >