shift的数据下载矩阵生成之路

要处理四个矩阵文件，分别为GSE109555 GSE136447 GSE36552 GSE66507

从ncbi网站上下载原始.sar文件，以及转化压缩成为.fastq.gz文件

在后台处理利用nohup&或者screen

cat GSE109555.txt | while read i;

/data/biosoft/software/sratoolkit.2.9.6-1-ubuntu64/bin/prefetch -X 9999999999999999999 -O `pwd` $i && echo "**${i}.sra done**";

fastq-dump --split-3 ${i}.sra;

gzip ${i}.fastq

done

————————————生成count矩阵，要求比对率为60%以上————————————

GSE66507，是双端测序文件，长度为50bp

比对过程中出现问题

多是由于数据没有下载完整造成的，此时需要重新下载

先是利用组里的index进行比对；结果还可以？大致有2/3比对率是大于百分之50的

GSE36552，是单端测序文件，长度为100bp

准确率较好，基本可以稳定在55%以上

但是对于

GSE109555，是双端测序文件，长度为150bp

以及

GSE136447，是双端测序文件，长度为150bp

比对率较差，都只有15%左右

更换比对文件Homo_sapiens.GRCh38.99.gtf;hg38_ercc.fa

STAR建立index文件

STAR --runMode genomeGenerate --genomeDir /data/shift/other/index_STAR --runThreadN 20 --genomeFastaFiles hg38_ercc.fa --sjdbGTFfile Homo_sapiens.GRCh38.99.2.gtf

出现问题：.fa和gtf文件里的chr1 和1之间有差异，求助解决

Solution: check the formatting of the GTF file. Most likely cause is the difference in chromosome naming between GTF and FASTA fill

确认文件中的开头

grep -o -E "^\w+([-+.]\w+)*" Homo_sapiens.GRCh38.99.gtf|sort|uniq

将.gtf文件中加chr

#! /user/bin/bash echo "usage: change_vcf_name.sh a | r input_vcf_file output_vcf_file" echo "a (add) : 添加chr字符到vcf文件中" echo "r (remove): 删除vcf文件中的chr字符"if [ $1 = "a" ]then awk '{ if($0 !~ /^#/) print "chr"$0; else if(match($0,/(##contig= $3elif [ $1 = "r" ]then sed 's/##contig= $3elseecho "输入的第一个参数不合适，请输入a或r。"fi

继续建立index文件，而后进行比对：本次比对对数据进行了清洗

#!/bin/bash

cat /data/shift/other/GSE36552/GSE36552.txt|while read i;

fastp -i /data/shift/other/GSE36552/GSE36552_data/${i}.fastq.gz -o /data/shift/other/GSE36552/GSE36552_data_clean/${i}_clean.fastq.gz;

STAR --runThreadN 12 --genomeDir /data/shift/other/index_STAR --readFilesCommand zcat --readFilesIn /data/shift/other/GSE36552/GSE36552_data_clean/${i}_clean.fastq.gz --sjdbGTFfile /data/shift/other/index_STAR/Homo_sapiens.GRCh38.99.2.gtf --outFileNamePrefix /data/shift/other/GSE36552/GSE36552_count.38.99/${i}_ --outSAMtype BAM SortedByCoordinate;

featureCounts -p -t exon -g gene_id -a /data/shift/other/index_STAR/Homo_sapiens.GRCh38.99.2.gtf -o /data/shift/other/GSE36552/GSE36552_count.38.99/${i}_counts.txt /data/shift/other/GSE36552/GSE36552_count.38.99/${i}_Aligned.sortedByCoord.out.bam;

done

但是比对结果并不尽如人意，比对成功率略低于组里文件：思考建立STAR index时是否应该根据不同的序列长度增加参数 --sjdbOverhang reads长度的最大值减1，默认是100

提出解决办法：随机选取一对文件fastq从中取出50bp(两个文件对称选取)的长度，与人进行比对

先对GSE136447尝试

利用软件seqkit

从Download - SeqKit - Ultrafast FASTA/Q kit (shenwei.me)下载软件

操作方法指路：seqkit：序列梳理神器-统计、格式转换、长度筛选、质量值转换、翻译、反向互补、抽样、去重、滑窗、拆分等30项全能..._刘永鑫Adam的博客-CSDN博客

报错：seqkit无法处理.gz文件

报错：更换剪切长度50，100都报错