测序数据比对到参考基因组

1.准备转录组数据

转录组数据可以从数据库中下载，也可以使用公司返回的测序数据，在服务器上下载完成后将原始文件名改为自己熟悉的格式，方便后期分析的进行。

2.准备样本信息表

样本信息表主要包括三部分，第一部分样品分组，第二部分：样本名称，第三部分：测序数据存在的绝对路径，例如以下的示例：

CB_S1   CB_S1_CF1   /path/to/your/workplace/CB_S1_CF1_1.fq.gz       /path/to/your/workplace/CB_S1_CF1_2.fq.gz
#第一列，样品分组。第二列，样本名称。第三列，正向测序文件存在的绝对位置。第四列，反向测序文件存在的绝对位置。

使用awk命令生成sample.txt文件

先将所有测序文件的文件名放在一个文件中，再替换掉不需要的部分，以下是相应代码。

ls *.gz > name.lst
#将所有以.gz结尾的文件的文件名写入一个名为name.lst的文件中
sed -i 's/.fq.gz//' name.lst
#将name.lst文件中所有.fq.gz替换成空格，sed -i 是原位替换，会直接更改文件中的内容
awk -F '_' '{print $1"_"$2"\t"$1"_"$2"_"$3"\t/path/to/your/workplace/"$1"_"$2"_"$3"_"$4".fq.gz""\t/path/to/your/workplace/"$1"_"$2"_"$3"_""2.fq.gz"}' name.lst
#使用awk命令批量生成代码，-F '_'表示以_为分隔符，\t为制表符，即大空格，其他部分同awk命令。

3. hisat2构建参考基因组

在正式比对前还需要构建参考基因组，所使用的软件是hisat2,基本命令为hisat2 -build,基础命令为：

hisat2-build /path/to/the/genome.fasta /path/to/your/output/genome 1>hisat2-build.log 2>&1

以上各代码分别为：
/path/to/the/genome.fasta：参考基因组所处位置;
/path/to/your/output/genome：输出文件所存储位置及所使用的前缀；
1>hisat2-build.log 2>&1：将标准输出流与错误输出流同时输入到hisat2_build.log这个文件中。
step1的结果文件会存放在参考基因组所在的文件夹

4.比对（参考基因组与测序结果生成sam文件）

记在前面的话：如果生成的比对日志文件比对率在70%以下，可能就需要调整命令中的相关参数

4.1 使用hisat2比对的基本命令

hisat2 --new-summary -p 6 -x /path/to/genome -U /path/to/the/seqdata -S outname.sam --rna-strandness R 1>name.log 2>&1

hisat2 --new-summary:命令的开头
-p 6：软件运行的线程数，一般为2-6
-x /path/to/genome：参考基因组所处的文件夹及构建的hisat2索引文件开头
-U /path/to/the/seqdata：测序文件所处的文件夹
**-S outname.sam **：比对结果的输出文件名
--rna-strandness：测序的策略，一般为非链特异性测序，所以此步可省略
R 1>name.log 2>&1：指定日志文件的输出的文件名

4.2 使用awk命令批量生成命令

要点：灵活使用常量与变量

awk '{print "hisat2 --new-summary -p 6 -x /path/to/genome -U "$3" -S "$2".sam --rna-strandness R 1>"$2".log 2>&1 &"}' sample.txt > step2.hisat_run.sh

以上命令是针对单末端测序项目的，单末端测序只产生一个测序文件，故直接使用-U后接$3即可，但是如果是双末端测序，则会产生两个测序文件，此时需要指定两个输入文件，命令有所变化，如下所示：

awk '{print "hisat2 --new-summary -p 6 -x /path/to/genome -1 "$3" -2 " $4" -S  "$2".sam --rna-strandness R 1>"$2".log 2>&1 &"}' sample.txt > step2.hisat_run.sh

4.3 samtools压缩文件生成bam文件

awk '{print "samtools sort -o "$2".bam "$2".sam"}' sample.txt > step3_sam2bam.sh

使用awk命令批量生成命令，使用的软件是samtools,生成的命令保存到step3_sam2bam.sh脚本中

4.4 samtools构建索引

awk '{print "samtools index "$2".bam"}' sample.txt > step4.bamindex.sh

使用awk命令批量生成命令，使用的软件是samtools,生成的命令保存到step4.bamindex.sh脚本中

5.IGV对测序比对结果可视化

5.1可视化所需数据

（1）基因组序列：参考基因组（genome.fasta；gene.gtf）
（2）压缩后的bam文件和对应的bai文件

5.2 IGV的使用

5.2.1创建基因组

Genomescreate a genome fileunique identifier(文件名)descriptive name(描述名称)fasta file(参考基因组)gene file(gtf文件)OK

5.2.1导入测序文件

fileload from filebam文件